---
title: 'Vicarious Contact PSRM Replication: Analysis'
author: "Lee-Or Ankori-Karlinsky, Rob Blair, Jessica Gottlieb, Samantha Moore-Berg"
output: 
  pdf_document:
    toc: yes
    toc_depth: 3
  html_document:
    toc: yes
    toc_depth: '3'
    df_print: paged
editor_options:
  chunk_output_type: console
---


```{r setup, include=FALSE}
#set Markdown Root directory / working directory 


####Set Appropriate Working Directory####
knitr::opts_knit$set(root.dir= "~/Brown Dropbox/Lee-Or Ankori-Karlinsky/Braver Angels Documentary/PSRM_ReplicationPackage/ContentContact_ReplicationPackage",echo = TRUE) 
setwd("~/Brown Dropbox/Lee-Or Ankori-Karlinsky/Braver Angels Documentary/PSRM_ReplicationPackage/ContentContact_ReplicationPackage") 

library(tidyverse)
library(haven)
library(dotwhisker)
library(ggpubr)
library(margins)
library(lmtest)
library(sandwich)
library(psych)
library(performance)
library(corrr)
library(ipw)
library(WeightIt)
library(coefplot)
library(usdata)
library(readxl)
library(kableExtra)
library(magrittr)
library(mice)
library(rlang)
library(estimatr)
library(texreg)
library(marginaleffects)
library(janitor)
library(fastDummies)
library(xtable)
library(patchwork)
library(hms)


#load cleaned data
final_drops <- read_csv("Data/Cleaned/final_drops.csv")


```


# Create Analysis + Plotting Functions: ATEs

```{r plotting functions}

####One Model ####

plot_generate_1 <- function(model_x, name_x, plot_title, xlim) {
  
  #create tidy df with robust SEs
  new_model_x <- coeftest(model_x, vcov=vcovHC, type="HC1") #robust SEs for model x
  name_x_df <- tidy(new_model_x) #create df
  name_x_df <- name_x_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_x}}) #remove intercept
  models <- name_x_df

  #plot
  plot <- dwplot(models,
        vline = geom_vline(
         xintercept = 0, 
         linetype= "dashed", 
         color = "black"),
       dot_args = list(size=5, shape=c(15)),
       whisker_args = list(size=0.7),
dodge_size = 0.8
) %>% 
relabel_predictors(treat_long_collapse2 = "Intention-To-Treat Effect 
                   of Vicarious Contact")+
  ggtitle(plot_title) +
  xlab("Coefficient Estimate with 95% Confidence Intervals")+
  xlim(xlim)+
  scale_colour_manual(
        name = "Legend",
        labels = list(models$model[1]),
        values = list("black")
        )+
  scale_shape_manual(
    name= "Legend",
    labels = list(models$model[1]),
    values=list(15)
  )+
  guides(color = guide_legend(override.aes=list(shape=c(15))))+
  theme_bw()+
  theme(plot.title = element_text(size=14, face="bold"), 
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size=12),
        legend.title = element_text(size=11, face="bold"),
        legend.text = element_text(size=10))

  #reduce to three decimals
  coef_1 <- format(round(c(models$estimate[1]), 3), nsmall=3)
  
  #add coefficient value
  plot + annotate(x=models$estimate[1], y=0.9, label=coef_1,geom="label", color="black")
  }


#Short videos
plot_generate_1_short <- function(model_x, name_x, plot_title, xlim) {
  
  #create tidy df with robust SEs
  new_model_x <- coeftest(model_x, vcov=vcovHC, type="HC1") #robust SEs for model x
  name_x_df <- tidy(new_model_x) #create df
  name_x_df <- name_x_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_x}}) #remove intercept
  models <- name_x_df

  #plot
  plot <- dwplot(models,
        vline = geom_vline(
         xintercept = 0, 
         linetype= "dashed", 
         color = "black"),
       dot_args = list(size=5, shape=c(15)),
       whisker_args = list(size=0.7),
dodge_size = 0.8
) %>% 
relabel_predictors(treat_short_collapse = "Short Videos vs.
                          Empty Control")+
  ggtitle(plot_title) +
  xlab("Coefficient Estimate with 95% Confidence Intervals")+
  xlim(xlim)+
  scale_colour_manual(
        name = "Legend",
        labels = list(models$model[1]),
        values = list("black")
        )+
  scale_shape_manual(
    name= "Legend",
    labels = list(models$model[1]),
    values=list(15)
  )+
  guides(color = guide_legend(override.aes=list(shape=c(15))))+
  theme_bw()+
  theme(plot.title = element_text(size=14, face="bold"), 
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size=12),
        legend.title = element_text(size=11, face="bold"),
        legend.text = element_text(size=10))

  #reduce to three decimals
  coef_1 <- format(round(c(models$estimate[1]), 3), nsmall=3)
  
  #add coefficient value
  plot + annotate(x=models$estimate[1], y=0.9, label=coef_1,geom="label", color="black")
  }



####Two Models####
plot_generate_2 <- function(model_x, model_y, name_x, name_y, plot_title, xlim) {
  
  #convert models into tidy dataframes with Robust SEs
  new_model_x <- coeftest(model_x, vcov=vcovHC, type="HC1") #robust SEs for model x
  new_model_y <- coeftest(model_y, vcov=vcovHC, type="HC1") #robust SEs for model y
  name_x_df <- tidy(new_model_x) #create df
  name_x_df <- name_x_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_x}}) #remove intercept
  name_y_df <- tidy(new_model_y) #create df
  name_y_df <- name_y_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_y}}) #remove intercept
  models <- rbind(name_x_df, name_y_df)

#plot
  plot <- dwplot(models,
        vline = geom_vline(
         xintercept = 0, 
         linetype= "dashed", 
         color = "black"),
       dot_args = list(size=5, shape=c(15,16)),
       whisker_args = list(size=0.7),
dodge_size = 0.8
) %>% 
relabel_predictors(treat_long_collapse2 = "Intention-To-Treat Effect 
                   of Vicarious Contact")+
  ggtitle(plot_title) +
  xlab("Coefficient Estimate with 95% Confidence Intervals")+
  xlim(xlim)+
  scale_colour_manual(
        name = "Legend",
        labels = list(models$model[2], models$model[1]),
        values = list("black", "black")
        )+
  scale_shape_manual(
    name= "Legend",
    labels = list(models$model[2], models$model[1]),
    values=list(15,16)
  )+
  guides(color = guide_legend(override.aes=list(shape=c(15:16))))+
  theme_bw()+
  theme(plot.title = element_text(size=14, face="bold"), 
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size=12),
        legend.title = element_text(size=11, face="bold"),
        legend.text = element_text(size=10))

  coef_1 <- format(round(c(models$estimate[1]), 3), nsmall=3)
  coef_2 <- format(round(c(models$estimate[2]), 3), nsmall=3)
  
  plot + annotate(x=models$estimate[2], y=1.3, label=coef_2,geom="label", color="black") +
  annotate(x=models$estimate[1], y=0.7, label=coef_1,geom="label", color="black")
  }

####Two Models No LABEL####
plot_generate_2_nolabel <- function(model_x, model_y, name_x, name_y, plot_title, xlim) {
  
  #convert models into tidy dataframes with Robust SEs
  new_model_x <- coeftest(model_x, vcov=vcovHC, type="HC1") #robust SEs for model x
  new_model_y <- coeftest(model_y, vcov=vcovHC, type="HC1") #robust SEs for model y
  name_x_df <- tidy(new_model_x) #create df
  name_x_df <- name_x_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_x}}) #remove intercept
  name_y_df <- tidy(new_model_y) #create df
  name_y_df <- name_y_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_y}}) #remove intercept
  models <- rbind(name_x_df, name_y_df)

#plot
  plot <- dwplot(models,
        vline = geom_vline(
         xintercept = 0, 
         linetype= "dashed", 
         color = "black"),
       dot_args = list(size=5, shape=c(15,16)),
       whisker_args = list(size=0.7),
dodge_size = 0.8
) %>% 
relabel_predictors(treat_long_collapse2 = "")+
  ggtitle(plot_title) +
  xlab("Coefficient Estimate with 95% Confidence Intervals")+
  xlim(xlim)+
  scale_colour_manual(
        name = "Legend",
        labels = list(models$model[2], models$model[1]),
        values = list("black", "black")
        )+
  scale_shape_manual(
    name= "Legend",
    labels = list(models$model[2], models$model[1]),
    values=list(15,16)
  )+
  guides(color = guide_legend(override.aes=list(shape=c(15:16))))+
  theme_bw()+
  theme(plot.title = element_text(size=14, face="bold"), 
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size=12),
        legend.title = element_text(size=11, face="bold"),
        legend.text = element_text(size=10))

  coef_1 <- format(round(c(models$estimate[1]), 3), nsmall=3)
  coef_2 <- format(round(c(models$estimate[2]), 3), nsmall=3)
  
  plot + annotate(x=models$estimate[2], y=1.3, label=coef_2,geom="label", color="black") +
  annotate(x=models$estimate[1], y=0.7, label=coef_1,geom="label", color="black")
  }




plot_generate_2_w3 <- function(model_x, model_y, name_x, name_y, plot_title, xlim) {
  
  #convert models into tidy dataframes with Robust SEs
  new_model_x <- coeftest(model_x, vcov=vcovHC, type="HC1") #robust SEs for model x
  new_model_y <- coeftest(model_y, vcov=vcovHC, type="HC1") #robust SEs for model y
  name_x_df <- tidy(new_model_x) #create df
  name_x_df <- name_x_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_x}}) #remove intercept
  name_y_df <- tidy(new_model_y) #create df
  name_y_df <- name_y_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_y}}) #remove intercept
  models <- rbind(name_x_df, name_y_df)

#plot
  plot <- dwplot(models,
        vline = geom_vline(
         xintercept = 0, 
         linetype= "dashed", 
         color = "black"),
       dot_args = list(size=5, shape=c(15,16)),
       whisker_args = list(size=0.7),
dodge_size = 0.8
) %>% 
relabel_predictors(treat_long_collapse2 = "Durability of Intention-To-Treat Effect 
                   of Vicarious Contact")+
  ggtitle(plot_title) +
  xlab("Coefficient Estimate with 95% Confidence Intervals")+
  xlim(xlim)+
  scale_colour_manual(
        name = "Legend",
        labels = list(models$model[2], models$model[1]),
        values = list("black", "black")
        )+
  scale_shape_manual(
    name= "Legend",
    labels = list(models$model[2], models$model[1]),
    values=list(15,16)
  )+
  guides(color = guide_legend(override.aes=list(shape=c(15:16))))+
  theme_bw()+
  theme(plot.title = element_text(size=14, face="bold"), 
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size=12),
        legend.title = element_text(size=11, face="bold"),
        legend.text = element_text(size=10))

  coef_1 <- format(round(c(models$estimate[1]), 3), nsmall=3)
  coef_2 <- format(round(c(models$estimate[2]), 3), nsmall=3)
  
  plot + annotate(x=models$estimate[2], y=1.3, label=coef_2,geom="label", color="black") +
  annotate(x=models$estimate[1], y=0.7, label=coef_1,geom="label", color="black")
  }


#short films
plot_generate_2_short <- function(model_x, model_y, name_x, name_y, plot_title, xlim) {
  
  #convert models into tidy dataframes with Robust SEs
  new_model_x <- coeftest(model_x, vcov=vcovHC, type="HC1") #robust SEs for model x
  new_model_y <- coeftest(model_y, vcov=vcovHC, type="HC1") #robust SEs for model y
  name_x_df <- tidy(new_model_x) #create df
  name_x_df <- name_x_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_x}}) #remove intercept
  name_y_df <- tidy(new_model_y) #create df
  name_y_df <- name_y_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_y}}) #remove intercept
  models <- rbind(name_x_df, name_y_df)

#plot
  plot <- dwplot(models,
        vline = geom_vline(
         xintercept = 0, 
         linetype= "dashed", 
         color = "black"),
       dot_args = list(size=5, shape=c(15,16)),
       whisker_args = list(size=0.7),
dodge_size = 0.8
) %>% 
relabel_predictors(treat_short_collapse = "Short Videos vs.
                          Empty Control")+
  ggtitle(plot_title) +
  xlab("Coefficient Estimate with 95% Confidence Intervals")+
  xlim(xlim)+
  scale_colour_manual(
        name = "Legend",
        labels = list(models$model[2], models$model[1]),
        values = list("black", "black")
        )+
  scale_shape_manual(
    name= "Legend",
    labels = list(models$model[2], models$model[1]),
    values=list(15,16)
  )+
  guides(color = guide_legend(override.aes=list(shape=c(15:16))))+
  theme_bw()+
  theme(plot.title = element_text(size=14, face="bold"), 
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size=12),
        legend.title = element_text(size=11, face="bold"),
        legend.text = element_text(size=10))

  coef_1 <- format(round(c(models$estimate[1]), 3), nsmall=3)
  coef_2 <- format(round(c(models$estimate[2]), 3), nsmall=3)
  
  plot + annotate(x=models$estimate[2], y=1.3, label=coef_2,geom="label", color="black") +
  annotate(x=models$estimate[1], y=0.7, label=coef_1,geom="label", color="black")
  }


####Three Models####
plot_generate_3 <- function(model_x, model_y, model_z, name_x, name_y, name_z, plot_title, xlim) {
  
  #convert models into tidy dataframes with Robust SEs
  new_model_x <- coeftest(model_x, vcov=vcovHC, type="HC1") #robust SEs for model x
  new_model_y <- coeftest(model_y, vcov=vcovHC, type="HC1") #robust SEs for model y
  new_model_z <- coeftest(model_z, vcov=vcovHC, type="HC1") #robust SEs for model z
  name_x_df <- tidy(new_model_x) #create df
  name_x_df <- name_x_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_x}}) #remove intercept
  name_y_df <- tidy(new_model_y) #create df
  name_y_df <- name_y_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_y}}) #remove intercept
  name_z_df <- tidy(new_model_z) #create df
  name_z_df <- name_z_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_z}})#remove intercept
  models <- rbind(name_x_df, name_y_df, name_z_df)

#plot
  plot <- dwplot(models,
        vline = geom_vline(
         xintercept = 0, 
         linetype= "dashed", 
         color = "black"),
       dot_args = list(size=5, shape=c(15,16,17)),
       whisker_args = list(size=0.7),
dodge_size = 0.8
) %>% 
relabel_predictors(treat_long_collapse2 = "Intention-To-Treat Effect 
                   of Vicarious Contact")+
  ggtitle(plot_title) +
  xlab("Coefficient Estimate with 95% Confidence Intervals")+
  xlim(xlim)+
  scale_colour_manual(
        name = "Legend",
        labels = list(models$model[3], models$model[2], models$model[1]),
        values = list("black", "black", "black")
        )+
  scale_shape_manual(
    name= "Legend",
    labels = list(models$model[3],models$model[2], models$model[1]),
    values=list(15,16, 17)
  )+
  guides(color = guide_legend(override.aes=list(shape=c(15:17))))+
  theme_bw()+
  theme(plot.title = element_text(size=14, face="bold"), 
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size=12),
        legend.title = element_text(size=11, face="bold"),
        legend.text = element_text(size=10))

  coef_1 <- format(round(c(models$estimate[1]), 3), nsmall=3)
  coef_2 <- format(round(c(models$estimate[2]), 3), nsmall=3)
  coef_3 <- format(round(c(models$estimate[3]), 3), nsmall=3)
  
  
  plot + annotate(x=models$estimate[3], y=1.35, label=coef_3,geom="label", color="black")+
    annotate(x=models$estimate[2], y=0.92, label=coef_2,geom="label", color="black") +
  annotate(x=models$estimate[1], y=0.65, label=coef_1,geom="label", color="black")
  }

####Three Models NO LABEL ####
plot_generate_3_nolabel <- function(model_x, model_y, model_z, name_x, name_y, name_z, plot_title, xlim) {
  
  #convert models into tidy dataframes with Robust SEs
  new_model_x <- coeftest(model_x, vcov=vcovHC, type="HC1") #robust SEs for model x
  new_model_y <- coeftest(model_y, vcov=vcovHC, type="HC1") #robust SEs for model y
  new_model_z <- coeftest(model_z, vcov=vcovHC, type="HC1") #robust SEs for model z
  name_x_df <- tidy(new_model_x) #create df
  name_x_df <- name_x_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_x}}) #remove intercept
  name_y_df <- tidy(new_model_y) #create df
  name_y_df <- name_y_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_y}}) #remove intercept
  name_z_df <- tidy(new_model_z) #create df
  name_z_df <- name_z_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_z}})#remove intercept
  models <- rbind(name_x_df, name_y_df, name_z_df)

#plot
  plot <- dwplot(models,
        vline = geom_vline(
         xintercept = 0, 
         linetype= "dashed", 
         color = "black"),
       dot_args = list(size=5, shape=c(15,16,17)),
       whisker_args = list(size=0.7),
dodge_size = 0.8
) %>% 
relabel_predictors(treat_long_collapse2 = "")+
  ggtitle(plot_title) +
  xlab("Coefficient Estimate with 95% Confidence Intervals")+
  xlim(xlim)+
  scale_colour_manual(
        name = "Legend",
        labels = list(models$model[3], models$model[2], models$model[1]),
        values = list("black", "black", "black")
        )+
  scale_shape_manual(
    name= "Legend",
    labels = list(models$model[3],models$model[2], models$model[1]),
    values=list(15,16, 17)
  )+
  guides(color = guide_legend(override.aes=list(shape=c(15:17))))+
  theme_bw()+
  theme(plot.title = element_text(size=14, face="bold"), 
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size=12),
        legend.title = element_text(size=11, face="bold"),
        legend.text = element_text(size=10))

  coef_1 <- format(round(c(models$estimate[1]), 3), nsmall=3)
  coef_2 <- format(round(c(models$estimate[2]), 3), nsmall=3)
  coef_3 <- format(round(c(models$estimate[3]), 3), nsmall=3)
  
  
  plot + annotate(x=models$estimate[3], y=1.35, label=coef_3,geom="label", color="black")+
    annotate(x=models$estimate[2], y=0.92, label=coef_2,geom="label", color="black") +
  annotate(x=models$estimate[1], y=0.65, label=coef_1,geom="label", color="black")
  }




####Four Models####
plot_generate_4 <- function(model_x, model_y, model_z, model_t, name_x, name_y, name_z, name_t, plot_title, xlim) {
  
  #convert models into tidy dataframes with Robust SEs
  new_model_x <- coeftest(model_x, vcov=vcovHC, type="HC1") #robust SEs for model x
  new_model_y <- coeftest(model_y, vcov=vcovHC, type="HC1") #robust SEs for model y
  new_model_z <- coeftest(model_z, vcov=vcovHC, type="HC1") #robust SEs for model z
  new_model_t <- coeftest(model_t, vcov=vcovHC, type="HC1") #robust SEs for model t
  name_x_df <- tidy(new_model_x) #create df
  name_x_df <- name_x_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_x}}) #remove intercept
  name_y_df <- tidy(new_model_y) #create df
  name_y_df <- name_y_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_y}}) #remove intercept
  name_z_df <- tidy(new_model_z) #create df
  name_z_df <- name_z_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_z}})#remove intercept
  name_t_df <- tidy(new_model_t) #create df
  name_t_df <- name_t_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_t}})#remove intercept
  models <- rbind(name_x_df, name_y_df, name_z_df, name_t_df)

#plot
  plot <- dwplot(models,
        vline = geom_vline(
         xintercept = 0, 
         linetype= "dashed", 
         color = "black"),
       dot_args = list(size=5, shape=c(15,16,17,18)),
       whisker_args = list(size=0.7),
dodge_size = 1
) %>% 
relabel_predictors(treat_long_collapse2 = "Intention-To-Treat Effect 
                   of Vicarious Contact")+
  ggtitle(plot_title) +
  xlab("Coefficient Estimate with 95% Confidence Intervals")+
  xlim(xlim)+
  scale_colour_manual(
        name = "Legend",
        labels = list(models$model[4], models$model[3], models$model[2], models$model[1]),
        values = list("black", "black", "black", "black")
        )+
  scale_shape_manual(
    name= "Legend",
    labels = list(models$model[4], models$model[3],models$model[2], models$model[1]),
    values=list(15,16, 17, 18)
  )+
  guides(color = guide_legend(override.aes=list(shape=c(15:18))))+
  theme_bw()+
  theme(plot.title = element_text(size=14, face="bold"), 
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size=12),
        legend.title = element_text(size=11, face="bold"),
        legend.text = element_text(size=10))

  coef_1 <- format(round(c(models$estimate[1]), 3), nsmall=3)
  coef_2 <- format(round(c(models$estimate[2]), 3), nsmall=3)
  coef_3 <- format(round(c(models$estimate[3]), 3), nsmall=3)
  coef_4 <- format(round(c(models$estimate[4]), 3), nsmall=3)
  
  
  plot + 
    annotate(x=models$estimate[4], y=1.45, label=coef_4,geom="label", color="black")+
    annotate(x=models$estimate[3], y=1.2, label=coef_3,geom="label", color="black")+
    annotate(x=models$estimate[2], y=0.8, label=coef_2,geom="label", color="black")+
    annotate(x=models$estimate[1], y=0.55, label=coef_1,geom="label", color="black")
  }

#### Four Models Alt####
plot_generate_4_2 <- function(model_x, model_y, model_z, model_t, name_x, name_y, name_z, name_t, plot_title, xlim) {
  
  #convert models into tidy dataframes with Robust SEs
  new_model_x <- coeftest(model_x, vcov=vcovHC, type="HC1") #robust SEs for model x
  new_model_y <- coeftest(model_y, vcov=vcovHC, type="HC1") #robust SEs for model y
  new_model_z <- coeftest(model_z, vcov=vcovHC, type="HC1") #robust SEs for model z
  new_model_t <- coeftest(model_t, vcov=vcovHC, type="HC1") #robust SEs for model t
  name_x_df <- tidy(new_model_x) #create df
  name_x_df <- name_x_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_x}}) #remove intercept
  name_y_df <- tidy(new_model_y) #create df
  name_y_df <- name_y_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_y}}) #remove intercept
  name_z_df <- tidy(new_model_z) #create df
  name_z_df <- name_z_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_z}})#remove intercept
  name_t_df <- tidy(new_model_t) #create df
  name_t_df <- name_t_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_t}})#remove intercept
  models <- rbind(name_x_df, name_y_df, name_z_df, name_t_df)

#plot
  plot <- dwplot(models,
        vline = geom_vline(
         xintercept = 0, 
         linetype= "dashed", 
         color = "black"),
       dot_args = list(size=5, shape=c(15,16,17,18)),
       whisker_args = list(size=0.7),
dodge_size = 1
) %>% 
relabel_predictors(treat_shortlong_collapse = "Intention-To-Treat Effect 
                   of Vicarious Contact")+
  ggtitle(plot_title) +
  xlab("Coefficient Estimate with 95% Confidence Intervals")+
  xlim(xlim)+
  scale_colour_manual(
        name = "Legend",
        labels = list(models$model[4], models$model[3], models$model[2], models$model[1]),
        values = list("black", "black", "black", "black")
        )+
  scale_shape_manual(
    name= "Legend",
    labels = list(models$model[4], models$model[3],models$model[2], models$model[1]),
    values=list(15,16, 17, 18)
  )+
  guides(color = guide_legend(override.aes=list(shape=c(15:18))))+
  theme_bw()+
  theme(plot.title = element_text(size=14, face="bold"), 
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size=12),
        legend.title = element_text(size=11, face="bold"),
        legend.text = element_text(size=10))

  coef_1 <- format(round(c(models$estimate[1]), 3), nsmall=3)
  coef_2 <- format(round(c(models$estimate[2]), 3), nsmall=3)
  coef_3 <- format(round(c(models$estimate[3]), 3), nsmall=3)
  coef_4 <- format(round(c(models$estimate[4]), 3), nsmall=3)
  
  
  plot + 
    annotate(x=models$estimate[4], y=1.45, label=coef_4,geom="label", color="black")+
    annotate(x=models$estimate[3], y=1.2, label=coef_3,geom="label", color="black")+
    annotate(x=models$estimate[2], y=0.8, label=coef_2,geom="label", color="black")+
    annotate(x=models$estimate[1], y=0.55, label=coef_1,geom="label", color="black")
  }


####Five Models####
plot_generate_5 <- function(model_x, model_y, model_z, model_t, model_w, name_x, name_y, name_z, name_t, name_w, plot_title, xlim) {
  
  #convert models into tidy dataframes with Robust SEs
  new_model_x <- coeftest(model_x, vcov=vcovHC, type="HC1") #robust SEs for model x
  new_model_y <- coeftest(model_y, vcov=vcovHC, type="HC1") #robust SEs for model y
  new_model_z <- coeftest(model_z, vcov=vcovHC, type="HC1") #robust SEs for model z
  new_model_t <- coeftest(model_t, vcov=vcovHC, type="HC1") #robust SEs for model t
  new_model_w <- coeftest(model_w, vcov=vcovHC, type="HC1") #robust SEs for model t
  name_x_df <- tidy(new_model_x) #create df
  name_x_df <- name_x_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_x}}) #remove intercept
  name_y_df <- tidy(new_model_y) #create df
  name_y_df <- name_y_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_y}}) #remove intercept
  name_z_df <- tidy(new_model_z) #create df
  name_z_df <- name_z_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_z}})#remove intercept
  name_t_df <- tidy(new_model_t) #create df
  name_t_df <- name_t_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_t}})#remove intercept
  name_w_df <- tidy(new_model_w) #create df
  name_w_df <- name_w_df %>% filter(term != "(Intercept)") %>% mutate(model = {{name_w}})#remove intercept
  models <- rbind(name_x_df, name_y_df, name_z_df, name_t_df, name_w_df)

#plot
  plot <- dwplot(models,
        vline = geom_vline(
         xintercept = 0, 
         linetype= "dashed", 
         color = "black"),
       dot_args = list(size=5, shape=c(15,16,17,18,19)),
       whisker_args = list(size=0.7),
dodge_size = 1
) %>% 
relabel_predictors(treat_long_collapse2 = "Intention-To-Treat Effect 
                   of Vicarious Contact")+
  ggtitle(plot_title) +
  xlab("Coefficient Estimate with 95% Confidence Intervals")+
  xlim(xlim)+
  scale_colour_manual(
        name = "Legend",
        labels = list(models$model[5], models$model[4], models$model[3], models$model[2], models$model[1]),
        values = list("black", "black", "black", "black", "black")
        )+
  scale_shape_manual(
    name= "Legend",
    labels = list(models$model[5], models$model[4], models$model[3],models$model[2], models$model[1]),
    values=list(15,16, 17, 18, 19)
  )+
  guides(color = guide_legend(override.aes=list(shape=c(15:19))))+
  theme_bw()+
  theme(plot.title = element_text(size=14, face="bold"), 
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size=12),
        legend.title = element_text(size=11, face="bold"),
        legend.text = element_text(size=10))

  coef_1 <- format(round(c(models$estimate[1]), 3), nsmall=3)
  coef_2 <- format(round(c(models$estimate[2]), 3), nsmall=3)
  coef_3 <- format(round(c(models$estimate[3]), 3), nsmall=3)
  coef_4 <- format(round(c(models$estimate[4]), 3), nsmall=3)
  coef_5 <- format(round(c(models$estimate[5]), 3), nsmall=3)
  
  
  plot + 
    annotate(x=models$estimate[5], y=1.48, label=coef_5,geom="label", color="black") +
    annotate(x=models$estimate[4], y=1.28, label=coef_4,geom="label", color="black")+
    annotate(x=models$estimate[3], y=1.08, label=coef_3,geom="label", color="black")+
    annotate(x=models$estimate[2], y=0.72, label=coef_2,geom="label", color="black")+
    annotate(x=models$estimate[1], y=0.52, label=coef_1,geom="label", color="black")
  }

####Six Models####
plot_generate_6 <- function(model_x, model_y, model_z, model_t, model_w, model_v, name_x, name_y, name_z, name_t, name_w, name_v, plot_title, xlim) {
  
  # Convert models into tidy dataframes with Robust SEs
  new_model_x <- coeftest(model_x, vcov=vcovHC, type="HC1") # robust SEs for model x
  new_model_y <- coeftest(model_y, vcov=vcovHC, type="HC1") # robust SEs for model y
  new_model_z <- coeftest(model_z, vcov=vcovHC, type="HC1") # robust SEs for model z
  new_model_t <- coeftest(model_t, vcov=vcovHC, type="HC1") # robust SEs for model t
  new_model_w <- coeftest(model_w, vcov=vcovHC, type="HC1") # robust SEs for model w
  new_model_v <- coeftest(model_v, vcov=vcovHC, type="HC1") # robust SEs for model v
  
  name_x_df <- tidy(new_model_x) %>% filter(term != "(Intercept)") %>% mutate(model = {{name_x}}) # remove intercept
  name_y_df <- tidy(new_model_y) %>% filter(term != "(Intercept)") %>% mutate(model = {{name_y}}) # remove intercept
  name_z_df <- tidy(new_model_z) %>% filter(term != "(Intercept)") %>% mutate(model = {{name_z}}) # remove intercept
  name_t_df <- tidy(new_model_t) %>% filter(term != "(Intercept)") %>% mutate(model = {{name_t}}) # remove intercept
  name_w_df <- tidy(new_model_w) %>% filter(term != "(Intercept)") %>% mutate(model = {{name_w}}) # remove intercept
  name_v_df <- tidy(new_model_v) %>% filter(term != "(Intercept)") %>% mutate(model = {{name_v}}) # remove intercept
  
  models <- rbind(name_x_df, name_y_df, name_z_df, name_t_df, name_w_df, name_v_df) # Combine models
  
  # Plot
  plot <- dwplot(models,
        vline = geom_vline(xintercept = 0, linetype = "dashed", color = "black"),
        dot_args = list(size = 5, shape = c(15, 16, 17, 18, 19, 20)),
        whisker_args = list(size = 0.7),
        dodge_size = 1
      ) %>%
      relabel_predictors(treat_long_collapse2 = "Intention-To-Treat Effect 
                         of Vicarious Contact") +
      ggtitle(plot_title) +
      xlab("Coefficient Estimate with 95% Confidence Intervals") +
      xlim(xlim) +
      scale_colour_manual(
        name = "Legend",
        labels = list(models$model[6], models$model[5], models$model[4], models$model[3], models$model[2], models$model[1]),
        values = list("black", "black", "black", "black", "black", "black")
      ) +
      scale_shape_manual(
        name = "Legend",
        labels = list(models$model[6], models$model[5], models$model[4], models$model[3], models$model[2], models$model[1]),
        values = list(15, 16, 17, 18, 19, 20)
      ) +
      guides(color = guide_legend(override.aes = list(shape = c(15:20)))) +
      theme_bw() +
      theme(
        plot.title = element_text(size = 14, face = "bold"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 11),
        axis.title.x = element_text(size = 12),
        legend.title = element_text(size = 11, face = "bold"),
        legend.text = element_text(size = 10)
      )
  
  coef_1 <- format(round(c(models$estimate[1]), 3), nsmall = 3)
  coef_2 <- format(round(c(models$estimate[2]), 3), nsmall = 3)
  coef_3 <- format(round(c(models$estimate[3]), 3), nsmall = 3)
  coef_4 <- format(round(c(models$estimate[4]), 3), nsmall = 3)
  coef_5 <- format(round(c(models$estimate[5]), 3), nsmall = 3)
  coef_6 <- format(round(c(models$estimate[6]), 3), nsmall = 3)
  
  plot + 
    annotate(x = models$estimate[6], y = 1.48, label = coef_6, geom = "label", color = "black") +
    annotate(x = models$estimate[5], y = 1.31, label = coef_5, geom = "label", color = "black") +
    annotate(x = models$estimate[4], y = 1.15, label = coef_4, geom = "label", color = "black") +
    annotate(x = models$estimate[3], y = 0.86, label = coef_3, geom = "label", color = "black") +
    annotate(x = models$estimate[2], y = 0.70, label = coef_2, geom = "label", color = "black") +
    annotate(x = models$estimate[1], y = 0.53, label = coef_1, geom = "label", color = "black")
}




####Seven Models####
plot_generate_7 <- function(models_list, names_list, plot_title, xlim) {
  
  models <- lapply(models_list, function(model) {
    coeftest(model, vcov = vcovHC, type = "HC1") %>%
      tidy() %>%
      filter(term != "(Intercept)")
  })
  
  names_df <- lapply(names_list, function(name) {
    data.frame(model = rep(name, nrow(models[[1]])))
  })
  
  models <- Map(cbind, models, names_df)
  
  models <- do.call(rbind, models)
  
  plot <- dwplot(models,
                 vline = geom_vline(
                   xintercept = 0,
                   linetype = "dashed",
                   color = "black"),
                 dot_args = list(size = 5),
                 whisker_args = list(size = 0.7),
                 dodge_size = 1
  ) %>%
    relabel_predictors(`treat_long_collapse2` = "Long Video vs. Placebo") +
    ggtitle(plot_title) +
    xlab("Coefficient Estimate with 95% Confidence Intervals") +
    xlim(xlim) +
    scale_colour_brewer(palette = "Dark2", name = "Model") +
    scale_shape_manual(values = rep(16, length(names_list))) +
    guides(color = rev(guide_legend(override.aes = list(shape = 16)))) +
    theme_minimal() +
    theme(
      plot.title = element_text(size = 14, face = "bold"),
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title.x = element_text(size = 12),
      legend.title = element_text(size = 11, face = "bold"),
      legend.text = element_text(size = 10)
    )
  
  plot
}



#### Seven Models Fix ####
plot_generate_7b <- function(models_list, names_list, plot_title, x_limits) {
  library(broom)
  library(sandwich)
  library(dotwhisker)
  library(dplyr)
  library(ggplot2)

  models <- lapply(models_list, function(model) {
    message("Tidying model")
    coeftest(model, vcov = vcovHC, type = "HC1") %>%
      tidy() %>%
      filter(term != "(Intercept)")
  })

  message("Tidied all models")

  # Make sure all models have the same number of terms
  term_counts <- sapply(models, nrow)
  if(length(unique(term_counts)) != 1) {
    stop("All models must have the same number of coefficients to bind correctly.")
  }

  names_df <- mapply(function(name, model_df) {
    data.frame(model = rep(name, nrow(model_df)))
  }, names_list, models, SIMPLIFY = FALSE)

  models <- Map(cbind, models, names_df)
  models_df <- do.call(rbind, models)

  message("Creating dwplot")

  p <- tryCatch({
    dwplot(models_df,
           vline = geom_vline(xintercept = 0, linetype = "dashed", color = "black"),
           dot_args = list(size = 5),
           whisker_args = list(size = 0.7),
           dodge_size = 1
    )
  }, error = function(e) {
    stop("dwplot failed: ", e$message)
  })

  message("After dwplot")
  print(class(p))

  p <- tryCatch({
    relabel_predictors(p, treat_long_collapse2 = "Long Video vs. Placebo")
  }, error = function(e) {
    stop("relabel_predictors failed: ", e$message)
  })

  message("After relabel_predictors")
  print(class(p))

  # Ensure x_limits is numeric
  if (!is.numeric(x_limits) || length(x_limits) != 2) {
    stop("x_limits must be a numeric vector of length 2.")
  }

  # Add aesthetics and return
  p <- p +
    ggtitle(plot_title) +
    xlab("Coefficient Estimate with 95% Confidence Intervals") +
    ggplot2::xlim(x_limits) +
    scale_colour_brewer(palette = "Dark2", name = "Model") +
    scale_shape_manual(values = rep(16, length(names_list))) +
    guides(color = guide_legend(override.aes = list(shape = 16))) +
    theme_minimal() +
    theme(
      plot.title = element_text(size = 14, face = "bold"),
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title.x = element_text(size = 12),
      legend.title = element_text(size = 11, face = "bold"),
      legend.text = element_text(size = 10)
    )

  message("Final plot built successfully.")
  return(p)
}






```

## Create Analysis + Plotting Functions: CATEs

```{r CATE plots}

####Create CATE DF####
get_cate <- function(model_x, name_1, name_2){
  model <- coeftest(model_x, vcov=vcovHC, type="HC1") #robust SEs for model_x
  model_df <- tidy(model) #create df
  model_df <- model_df %>% #remove intercept
    filter(term != "(Intercept)")
  
  #get estimate (Beta 1 + Beta 3)
  estimate = model_df$estimate[1] + model_df$estimate[3] 
  
  #get standard errors using vcov matrix
  vc <- vcovHC(model_x, type="HC1")
  varb1 <- vc[2, 2]
  varb3 <- vc[4, 4]
  covb1b3 <- vc[2, 4]
  se <- sqrt(varb1+varb3+2*covb1b3)
  
  #get degrees of freedom
  df_res <- model_x$df.residual
  
  #insert values df
  model_df <- model_df %>% 
    add_row(term={{name_2}}, 
            estimate = estimate, 
            std.error = se, 
            statistic = estimate/se, #t stat
            p.value = 2*pt(abs(statistic), df_res , lower.tail=F)
    ) %>% 
    slice(1,3,4) %>% #reorder variables
    mutate(term = c({{name_1}}, "Interaction", {{name_2}}) #rename terms
           )
}



#no interaction
get_cate_2 <- function(model_x, name_1, name_2){
  model <- coeftest(model_x, vcov=vcovHC, type="HC1") #robust SEs for model_x
  model_df <- tidy(model) #create df
  model_df <- model_df %>% #remove intercept
    filter(term != "(Intercept)")
  
  #get estimate (Beta 1 + Beta 3)
  estimate = model_df$estimate[1]+ model_df$estimate[3] 
  
  #get standard errors using vcov matrix
  vc <- vcovHC(model_x, type="HC1")
  varb1 <- vc[2, 2]
  varb3 <- vc[4, 4]
  covb1b3 <- vc[2, 4]
  se <- sqrt(varb1+varb3+2*covb1b3)
  
  #get degrees of freedom
  df_res <- model_x$df.residual
  
  #insert values df
  model_df <- model_df %>% 
    add_row(term={{name_2}}, 
            estimate = estimate, 
            std.error = se, 
            statistic = estimate/se, #t stat
            p.value = 2*pt(abs(statistic), df_res , lower.tail=F)
    ) %>% 
    slice(1,4) %>% #reorder variables
    mutate(term = c({{name_1}}, {{name_2}}) #rename terms
           )
}


get_cate_covars <- function(model_x, name_1, name_2, covariates = NULL) {
  model <- coeftest(model_x, vcov = vcovHC, type = "HC1") # robust SEs
  model_df <- tidy(model) # create df
  model_df <- model_df %>% 
  filter(!term %in% c("(Intercept)", covars))
  
  # Calculate estimate and SE for interaction term
  estimate <- model_df$estimate[1] + model_df$estimate[3]
  vc <- vcovHC(model_x, type = "HC1")
  varb1 <- vc[2, 2]
  varb3 <- vc[4, 4]
  covb1b3 <- vc[2, 4]
  se <- sqrt(varb1 + varb3 + 2 * covb1b3)
  df_res <- model_x$df.residual
  
  # Create the final data frame
  model_df <- model_df %>% 
    add_row(
      term = {{name_2}}, 
      estimate = estimate, 
      std.error = se, 
      statistic = estimate / se, 
      p.value = 2 * pt(abs(estimate / se), df_res, lower.tail = FALSE)
    ) %>% 
    slice(1, 3, 4) %>% 
    mutate(term = c({{name_1}}, "Interaction", {{name_2}}))
  return(model_df)
}


####Plot CATE####

plot_cate <- function(cate_df, plot_title, xlim){
  
  #name df
  df <- cate_df %>% 
    filter(term !="Interaction")
  
  #plot
  plot <- dwplot(df, 
         vline = geom_vline(
         xintercept = 0, 
         linetype= "dashed", 
         color = "black"),
       dot_args = list(size=5, color="black"),
       whisker_args = list(size=0.7, color="black"),
       dodge_size = 1) +
  ggtitle(plot_title)+
  xlab("Coefficient Estimate with 95% Confidence Intervals")+
  xlim(xlim)+
  theme_bw() + 
  theme(plot.title = element_text(size=14, face="bold"), 
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size=12),
        legend.title = element_text(size=11, face="bold"),
        legend.text = element_text(size=10))
  
  #create HTE label
  
  HTE_estimate <- format(round(c(cate_df$estimate[2]), 3), nsmall=3)
  HTE_pvalue <- format(round(c(cate_df$p.value[2]), 3), nsmall=3)
  HTE_label <- paste0("HTE: beta=", HTE_estimate, ", ", 
                      "p.value=", HTE_pvalue)
  
  plot + 
    annotate(geom = "label", label = HTE_label, x=0.3, y = 0.7, size = 4)
}

plot_cate2 <- function(cate_df, plot_title, xlim){
  
  # Extract estimate and p-value for HTE label
  HTE_estimate <- format(round(c(cate_df$estimate[2]), 3), nsmall=3)
  HTE_pvalue <- format(round(c(cate_df$p.value[2]), 3), nsmall=3)
  HTE_label <- paste0("HTE: beta=", HTE_estimate, ", p.value=", HTE_pvalue)
  
  # Name df and plot
  df <- cate_df %>% 
    filter(term != "Interaction")
  
  plot <- dwplot(df, 
         vline = geom_vline(
         xintercept = 0, 
         linetype= "dashed", 
         color = "black"),
       dot_args = list(size=5, color="black"),
       whisker_args = list(size=0.7, color="black"),
       dodge_size = 1) +
  ggtitle(plot_title)+
  xlab("Coefficient Estimate with 95% Confidence Intervals")+
  xlim(xlim)+
  theme_bw() + 
  theme(plot.title = element_text(size=14, face="bold"), 
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size=12),
        legend.title = element_text(size=11, face="bold"),
        legend.text = element_text(size=10))
  
  # Add HTE label to the plot
  plot <- plot + 
    annotate(geom = "label", label = HTE_label, x = 0.3, y = 0.7, size = 4)
  
  # Return the plot and the label (HTE_label) for later use
  return(list(plot = plot, label = HTE_label))
}



```


## Regression Tables Function

```{r create tables}

#Table functions take lm_robust models as input (robust standard errors come with the model) 

#### One Model####
get_table_1 <- function(model, table_title, model_name){
  texreg(
    model,
    custom.model.names = model_name,
    caption = table_title,
    caption.above = T,
    omit.coef = "(Intercept)",
    include.ci=F,
    fontsize = "large",
    float.pos = "H",
    threeparttable = T, 
    custom.note = "\\item %stars. \\\\ Robust standard errors are used."
  )
  }


#### Two Models ####
get_table_2 <- function(model1, model2, table_title, model1_name, model2_name){
  
  models <- list(model1, model2)
  
  texreg(
    models,
    custom.model.names = c(model1_name, model2_name),
    caption = table_title,
    caption.above = T,
    omit.coef = "(Intercept)",
    include.ci=F,
    fontsize = "large",
    float.pos = "H",
    threeparttable = T, 
    custom.note = "\\item %stars. \\\\ Robust standard errors are used."
  )
  }



#### Three Models ####
get_table_3 <- function(model1, model2, model3, table_title, model1_name, model2_name, model3_name){
  
  models <- list(model1, model2, model3)
  
  texreg(
    models,
    custom.model.names = c(model1_name, model2_name, model3_name),
    caption = table_title,
    caption.above = T,
    omit.coef = "(Intercept)",
    include.ci=F,
    fontsize = "large",
    float.pos = "H",
    threeparttable = T, 
    custom.note = "\\item %stars. \\\\ Robust standard errors are used."
  )
  }


####Four Models####
get_table_4 <- function(model1, model2, model3, model4, table_title, model1_name, model2_name, model3_name, model4_name){
  
  models <- list(model1, model2, model3, model4)
  
  texreg(
    models,
    custom.model.names = c(model1_name, model2_name, model3_name, model4_name),
    caption = table_title,
    caption.above = T,
    omit.coef = "(Intercept)",
    include.ci=F,
    fontsize = "large",
    float.pos = "H",
    threeparttable = T, 
    custom.note = "\\item %stars. \\\\ Robust standard errors are used."
  )
  }

```


## Model Summary Lines for Latex

```{r stat summary line}

#Regular 

latex_summary <- function(model) {
  # Compute robust standard errors (HC1)
  robust_se <- sqrt(diag(vcovHC(model, type = "HC1")))
  
  # Include the coefficient estimate
  coef_est <- coef(model)[-1]
  
  # Construct robust confidence intervals
  ci_str <- paste0(", 95\\% CI ", format(round(coef_est - 1.96 * robust_se, 3), nsmall = 3), " to ", format(round(coef_est + 1.96 * robust_se, 3), nsmall = 3), ",")
  
  # Calculate p-value using robust standard errors
  robust_p_value <- 2 * (1 - pnorm(abs(coef_est) / robust_se))
  p_value_str <- format(round(robust_p_value, 3), nsmall = 3)
  
  # Get the number of observations
  n_obs <- length(model$model[, 1])
  n_obs_str <- paste0("$N$ = ", n_obs)
  
  latex_line <- paste("$\beta$ =", format(round(coef_est, 3), nsmall = 3), ci_str, "$P$ =", p_value_str, ",", n_obs_str)
  
  return(latex_line)
}


#CATE
latex_summary_cate <- function(df, model) {
  output_lines <- df %>%
    rowwise() %>%
    mutate(
      term = paste("$\beta$ =", format(round(estimate, 3), nsmall = 3)),
      ci_str = paste0("95\\% CI", format(round(estimate - 1.96 * std.error, 3), nsmall = 3), " to ", format(round(estimate + 1.96 * std.error, 3), nsmall = 3)),
      p_value_str = paste0("$P$=", format(round(p.value,3), nsmall = 3))
    ) %>%
    select(term, ci_str, p_value_str) %>%
    summarise(latex_line = paste(term, ci_str, p_value_str, collapse = " "))

  # Extract the number of observations from the model object
  n_obs <- nobs(model)
  n_obs_str <- paste0("$N$=", n_obs)

  output_lines$latex_line <- paste(output_lines$latex_line, n_obs_str)

  return(output_lines$latex_line)
}



```


# Analysis 


## Figure 1: Vicarious intergroup contact reduced affective polarization among a nationally representative sample of Americans


```{r Figure 1}

####Regression####

#placebo vs long (with IPWs) 
H1.1 <- lm(aff_pol_idx~treat_long_collapse2, final_drops, weights=ipw_matched_w2_noimp)
#summary(H1.1)
#coeftest(H1.1, vcov = vcovHC, type="HC1") #Validated using Stata on 3/29/23 by RB


#placebo vs. long: outparty only (with IPWs)
H1.1_out <- lm(aff_pol_idx_outparty~treat_long_collapse2, final_drops, weights=ipw_matched_w2_noimp)
#summary(H1.1_out)
#coeftest(H1.1_out, vcov = vcovHC, type="HC1")

#summary statistics for in-text references
latex_summary(H1.1)
latex_summary(H1.1_out)

####Plot####
Main_AffPol_plot <- plot_generate_2(H1.1, H1.1_out, "Affective Polarization", "Outparty Only", "", c(-0.4, 0.1))


#save to Figures folder
ggsave("Figures/Figure1.png", plot = Main_AffPol_plot, width=10, height=7)


```



## Figure 2: Vicarious contact increases interest but not investment in depolarization activities

```{r Figure 2}

##### Regressions ####

#BA newsletter click 
H1.2a_click <- lm(BA_newsletter_clicked~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)
#summary(H1.2a_click)
#coeftest(H1.2a_click, vcov = vcovHC, type="HC1") #Validated using Stata on 3/29/23 by RB

latex_summary(H1.2a_click)

#donations 
H1.2b <- lm(donate_any~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)
#summary(H1.2b)
#coeftest(H1.2b, vcov = vcovHC, type="HC1") #Validated using Stata on 3/29/23 by RB

latex_summary(H1.2b)

#### Plot ####

#investment and interest plot 
behavior_plot <- plot_generate_2(H1.2a_click, H1.2b, "Interest in Depolarization", "Investment in Depolarization", "", c(-0.2, 0.2))

#show
print(behavior_plot)

#save to figures folder
ggsave("Figures/Figure2.png", plot = behavior_plot, width=10, height=7)

```


## Figure 3: Vicarious contact increased optimism and strengthened belief in the efficacy of dialogue

```{r Figure 3}

#### Regression ####

#belief in possibility of non-violent democratic change 
Y1.1 <- lm(nonviol_change_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)

#belief in efficacy of dialogue 
Y1.4 <- lm(dialogue_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)

#optimism about restoring civility
Y2.1 <- lm(optimism_civil_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)

#optimism about survival of democratic institutions 
Y2.2 <- lm(optimism_survive_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)

latex_summary(Y1.1)
latex_summary(Y1.4)
latex_summary(Y2.1)
latex_summary(Y2.2)

#### Plot ####

#New plot, created 11/26 for paper 
optimism_plot_PAP <- plot_generate_4(Y1.1, Y1.4, Y2.1, Y2.2, "Non-Violent Change Still Possible", "Belief in Efficacy of Dialogue","Restoring Civility", "Survival of Democratic Institutions","", c(-0.2, .7))

print(optimism_plot_PAP)

#save
ggsave("Figures/Figure3.png", plot = optimism_plot_PAP, width=10, height=7)

```

## Figure 4: Vicarious contact does not affect support for anti-democratic actions

```{r Figure 4}

####Regression####
H2.1 <- lm(antidem_idx~treat_long_collapse2, final_drops, weights=ipw_matched_w2_noimp)

latex_summary(H2.1)

Antidem_main_plot <- plot_generate_1(H2.1, "Anti-Democratic Index", "", c(-.3, .1))

#check
print(Antidem_main_plot)

#save plot
ggsave("Figures/Figure4.png", plot = Antidem_main_plot, width=10, height=7)


```

## Figure 5: Vicarious contact reduced affective polarization primarily among Democrats

```{r Figure 5}

####Regressions####
HTE1.1 <- lm(aff_pol_idx~treat_long_collapse2*partyID, final_drops, weights = ipw_matched_w2_noimp)

#get df for cate
HTE1.1_df <- get_cate(HTE1.1, "Democrats", "Republicans")

head(HTE1.1_df)

summary(HTE1.1)

#stats table
xtable::xtable(HTE1.1_df, "CATE Statistics: Affective Polarization")


#plot
AffPol_CATE_partyID <- plot_cate(HTE1.1_df, "", c(-0.4, 0.4))

#check
print(AffPol_CATE_partyID)

#save plot
ggsave("Figures/Figure5.png", plot = AffPol_CATE_partyID, width=10, height=7)


```

## Figure 6 

```{r Figure 6}

####Regression####
#placebo vs long (with IPWs) 
R3.1 <- lm(aff_pol_idx_w3~treat_long_collapse2, final_drops, weights=ipw_matched_w3_noimp)

#placebo vs. long: outparty only (with IPWs)
R3.1_out <- lm(aff_pol_idx_outparty_w3~treat_long_collapse2, final_drops, weights=ipw_matched_w3_noimp)

latex_summary(R3.1)
latex_summary(R3.1_out)


####Plot####
W3_AffPol_plot <- plot_generate_2_w3(R3.1, R3.1_out, "Affective Polarization", "Outparty Only", "", c(-0.4, 0.1))

#check
print(W3_AffPol_plot)

#save
ggsave("Figures/Figure6.png", plot = W3_AffPol_plot, width=10, height=7)

```


# Supplementary Information Tables and Figures 

## Descriptive statistics 


```{r demographics (Table SI.1)}

demos <- final_drops %>% 
  group_by(video_treatment) %>% 
  summarise(N = length(video_treatment),
            mean_age = mean(age), 
            pct_college = (sum(educ_dum)/N)*100, #4-year or more=1
            pct_female =  (sum(sex)/N)*100, #male is 0
            pct_republican = (sum(partyID)/N)*100,
            pct_child = sum(child/N, na.rm = T)*100, #child
            pct_white = sum(white/N, na.rm = T)*100,
            pct_christian = sum(christian/N, na.rm = T)*100
  )

demos_2 <- final_drops %>% 
  summarise(
            mean_age = mean(age), 
            pct_college = sum(educ_dum/2573)*100, #4-year or more=1
            pct_female =  sum(sex/2573)*100, #male is 0
            pct_republican = sum(partyID/2573)*100,
            pct_child = sum(child/2573, na.rm = T)*100, #child
            pct_white = sum(white/2573, na.rm = T)*100,
            pct_christian = sum(christian/2573, na.rm = T)*100
  ) %>% 
    mutate(
    video_treatment = "Total",
    N = 2573) %>% 
  select(video_treatment, N, mean_age, pct_college, pct_female, pct_republican, pct_child, pct_white, pct_christian)

demos <- rbind(demos, demos_2)

####Means and SDs####

#create new variables 
descriptives <- final_drops %>% 
  mutate(
    Republican = if_else(partyID == 1,1,0),
    Democrat = if_else(partyID == 0,1,0),
    Black = if_else(race_2 == 1, 1,0),
    Hispanic = if_else(race_2 == 2, 1,0),
    Asian = if_else(race_2 == 4, 1,0),
    Male = if_else(sex == 0, 1, 0),
    Female = if_else(sex==1,1,0)
    )

# List of variables
variable_list <- list(
  descriptives$Republican, 
  descriptives$Democrat, 
  descriptives$age, 
  descriptives$Male, 
  descriptives$Female, 
  descriptives$white, 
  descriptives$Black, 
  descriptives$Hispanic, 
  descriptives$Asian, 
  descriptives$educ_dum
)

# List of variable names
variable_names <- c(
  "Republican", "Democrat", "Age", "Male", "Female",
  "White", "Black", "Hispanic", "Asian", "College"
)

# Function to calculate descriptive statistics
calculate_stats <- function(variable_name, variable) {
  n <- ifelse(variable_name == "Age", sum(!is.na(variable)), sum(variable))
  mean_val <- mean(variable)
  sd_val <- sd(variable)
  return(c(n, format(mean_val, digits = 2), format(sd_val, digits = 2)))
}

# Calculate statistics for each variable
stats_list <- mapply(calculate_stats, variable_names, variable_list, SIMPLIFY = FALSE)

# Convert the list of statistics to a matrix and add column names
stats_matrix <- do.call(rbind, stats_list)
col_names <- c("N", "Mean", "SD")
rownames(stats_matrix) <- variable_names
colnames(stats_matrix) <- col_names

# Create a data frame for the results
stats_df <- as.data.frame(stats_matrix, stringsAsFactors = FALSE)



####Nationally representative (matched)####

d2 <- descriptives %>% 
  dplyr::filter(!is.na(weight))

# List of variables
variable_list <- list(
  d2$Republican, 
  d2$Democrat, 
  d2$age, 
  d2$Male, 
  d2$Female, 
  d2$white, 
  d2$Black, 
  d2$Hispanic, 
  d2$Asian, 
  d2$educ_dum
)


# Calculate statistics for each variable
stats_list2 <- mapply(calculate_stats, variable_names, variable_list, SIMPLIFY = FALSE)

# Convert the list of statistics to a matrix and add column names
stats_matrix2 <- do.call(rbind, stats_list2)
col_names <- c("N", "Mean", "SD")
rownames(stats_matrix2) <- variable_names
colnames(stats_matrix2) <- col_names

# Create a data frame for the results
stats_df2 <- as.data.frame(stats_matrix2, stringsAsFactors = FALSE)

#combine
descriptive_stats <- cbind(stats_df2, stats_df)

#Tables
xtable::xtable(descriptive_stats, caption = "Descriptive statistics")

```

```{r descriptive statistics outcomes (Table SI.2)}

# Variable name
# List of variables
variable_list_w2 <- list(
  descriptives$aff_pol_idx, 
  descriptives$aff_pol_idx_outparty, 
  descriptives$BA_newsletter_clicked, 
  descriptives$donate_all, 
  descriptives$antidem_idx, 
  descriptives$stereo_outparty_neg, 
  descriptives$stereo_outparty_pos, 
  descriptives$mass_abortion_diff,
  descriptives$mass_marriage_diff,
  descriptives$mass_leave_diff,
  descriptives$optimism_survive, 
  descriptives$optimism_pol,
  descriptives$optimism_civil,
  descriptives$dialogue
)


# List of variable names
DV_names <- c(
  "Affective Polarization", "Affective Polarization (Outparty Animus)", "BA Newsletter Clicks", "All Donations", "Anti-Democratic Attitudes",
  "Positive Outparty Stereotypes", "Negative Outparty Stereotypes", "Mass Perceptions: Abortion", "Mass Perceptions: Marriage", "Mass Perceptions: Parental Leave", "Optimism: survival of demcoratic institutions", "Optimism: overcoming polarization", "Optimism: restoring civility", "Dialogue as effective tool for change")

# Function to calculate descriptive statistics
calculate_DVs <- function(variable_name, variable) {
  n <- sum(!is.na(variable))
  mean_val <- mean(variable, na.rm = T)
  sd_val <- sd(variable, na.rm = T)
  return(c(n, format(mean_val, digits = 2), format(sd_val, digits = 2)))
}

# Calculate statistics for each variable
stats_list_3 <- mapply(calculate_DVs, DV_names, variable_list_w2, SIMPLIFY = FALSE)

# Convert the list of statistics to a matrix and add column names
stats_matrix_3 <- do.call(rbind, stats_list_3)
col_names <- c("N", "Mean", "SD")
rownames(stats_matrix_3) <- DV_names
colnames(stats_matrix_3) <- col_names

# Create a data frame for the results
DV_df_3 <- as.data.frame(stats_matrix_3, stringsAsFactors = FALSE)

DV_df_3 <- DV_df_3 %>% 
  select("Mean", "SD", "N")


#### Nationally representative sample####
# List of variables
variable_list_w2_2 <- list(
  d2$aff_pol_idx, 
  d2$aff_pol_idx_outparty, 
  d2$BA_newsletter_clicked, 
  d2$donate_all, 
  d2$antidem_idx, 
  d2$stereo_outparty_neg, 
  d2$stereo_outparty_pos, 
  d2$mass_abortion_diff,
  d2$mass_marriage_diff,
  d2$mass_leave_diff,
  d2$optimism_survive, 
  d2$optimism_pol,
  d2$optimism_civil,
  d2$dialogue
)


# Function to calculate descriptive statistics
calculate_DVs <- function(variable_name, variable) {
  n <- sum(!is.na(variable))
  mean_val <- mean(variable, na.rm = T)
  sd_val <- sd(variable, na.rm = T)
  return(c(n, format(mean_val, digits = 2), format(sd_val, digits = 2)))
}

# Calculate statistics for each variable
stats_list_4 <- mapply(calculate_DVs, DV_names, variable_list_w2_2, SIMPLIFY = FALSE)

# Convert the list of statistics to a matrix and add column names
stats_matrix_4 <- do.call(rbind, stats_list_4)
col_names <- c("N", "Mean", "SD")
rownames(stats_matrix_4) <- DV_names
colnames(stats_matrix_4) <- col_names

# Create a data frame for the results
DV_df_4 <- as.data.frame(stats_matrix_4, stringsAsFactors = FALSE)

DV_df_4 <- DV_df_4 %>% 
  select("Mean", "SD", "N")

#combine
DV_descriptives <- cbind(DV_df_4, DV_df_3)

xtable::xtable(DV_descriptives)


```

### Participants by condition

```{r participants by condition (Tables SI.3-4)}

#nationally representative (Table SI.3)
matched_sample <- final_drops %>% 
  dplyr::filter(!is.na(weight))

treatment2 <- matched_sample %>% 
  group_by(video_treatment) %>% 
  summarise(Randomized = length(video_treatment), 
            Wave2_Attriters = sum(matched_attriter_w2),
            N_Wave2 = Randomized-Wave2_Attriters,
            Wave3_Attriters = sum(matched_attriter_w3),
            N_Wave3 = Randomized-Wave3_Attriters)

xtable::xtable(treatment2, digits = 0)


#full sample (Table SI.4)
treatment <- final_drops %>% 
  group_by(video_treatment) %>% 
  summarise(Randomized = length(video_treatment), 
            Wave2_Attriters = sum(rand_attriter),
            N_Wave2 = Randomized-Wave2_Attriters,
            Wave3_Attriters = sum(rand_attriter_w3),
            N_Wave3 = Randomized-Wave3_Attriters)

xtable::xtable(treatment, digits = 0)


```

### Correlation Tables 

```{r correlation tables (Tables SI.5-6)}

#pre-registered index
index_items <- final_drops[, c("trust_diff_scale", "therm_diff_scale", "discomfort_outparty_friends_scale", "discomfort_outparty_marriage_scale", "discomfort_outparty_neighbors_scale", "threat_outparty_scale", "neg_partisan_scale")]


#change names of index items 
colnames(index_items) <- c("Trust", "FT", "Friends", "Marriage", "Neighbors", "Threat", "Negative_Par")


#confirm Cronbach's alpha is 0.74 (confirmed by LAK, 11/24/24)
alpha(index_items)

#outparty only 
index_items_out <- final_drops[, c("trust_outparty_scale", "therm_outparty_scale", "discomfort_outparty_friends_scale", "discomfort_outparty_marriage_scale", "discomfort_outparty_neighbors_scale", "threat_outparty_scale", "neg_partisan_scale")]

#change names of index items 
colnames(index_items_out) <- c("Trust (out)", "FT (out)", "Friends", "Marriage", "Neighbors", "Threat", "Negative_Par")

#confirm Cronbach's alpha is 0.77 (confirmed by LAK, 11/24/24; rounded up here)
alpha(index_items_out)


#regular index
aff_polidx_cor_matrix <- cor(index_items, use = "pairwise.complete.obs")  # Use pairwise deletion for missing values

#round
aff_polidx_cor_matrix <- round(aff_polidx_cor_matrix, 3)

print(aff_polidx_cor_matrix)

xtable::xtable(aff_polidx_cor_matrix)


#outparty index 
aff_polidx_out_cor_matrix <- cor(index_items_out, use = "pairwise.complete.obs")  # Use pairwise deletion for missing values

aff_polidx_out_cor_matrix <- round(aff_polidx_out_cor_matrix, 3)

print(aff_polidx_out_cor_matrix)

xtable::xtable(aff_polidx_out_cor_matrix)


```


## Balance tests 

```{r Balance tests (SI.7-8)}

long_placebo <- final_drops %>% 
  dplyr::filter(!is.na(treat_long_collapse2)) %>% 
  mutate(
    matched = if_else(!is.na(weight), 1, 0)
  )


full_model <- lm(treat_long_collapse2~partyID+age+white+Black+Asian+Hispanic+Female+educ_dum+ideology_yg, descriptives)

full_model_2 <- lm(treat_all_collapse~partyID+age+white+Black+Asian+Hispanic+Female+educ_dum+ideology_yg, descriptives)

full_model_3 <- lm(treat_long_collapse2~partyID+age+white+Black+Asian+Hispanic+Female+educ_dum+ideology_yg, d2)

full_model_4 <- lm(treat_all_collapse~partyID+age+white+Black+Asian+Hispanic+Female+educ_dum+ideology_yg, d2)


stargazer::stargazer(list(full_model_3, full_model_4), omit = "Constant", title = "Balance Test: Assignment to Treatment as Predicted by Covariates")

stargazer::stargazer(list(full_model, full_model_2), omit = "Constant", title = "Balance Test: Assignment to Treatment as Predicted by Covariates")



```

## Attrition tables 

```{r tables SI.9-11}

#### Table SI.9 ####

# Function to extract OLS estimates and format results
extract_ols_results <- function(model) {
  coef_summary <- summary(model)$coefficients
  estimate <- coef_summary[2, 1]
  p_value <- coef_summary[2, 4]
  
  # Add stars based on p-value
  stars <- case_when(
    p_value < 0.001 ~ "***",
    p_value < 0.01  ~ "**",
    p_value < 0.05  ~ "*",
    TRUE            ~ ""
  )
  
  # Return formatted string with estimate and stars
  sprintf("%.3f (p=%.3f)%s", estimate, p_value, stars)
}

# Run models for different conditions and waves (example)
conditions <- c("Full Film vs. Placebo", "Full Film vs. Placebo + Control", "Full Film vs. Short PMC", "Full Film vs. Short VC", "Full Film vs. Both Short Videos", "Short PMC vs. Empty Control", "Short VC vs. Empty Control", "Both Short Videos vs. Empty Control", "Short VC vs. Short PMC") # Add more comparisons here
models_w2 <- list(
  lm(matched_attriter_w2 ~ treat_long_collapse2, final_drops), #long vs. placebo
  lm(matched_attriter_w2 ~ treat_long_collapse1, final_drops), #long vs. placebo+control 
  lm(matched_attriter_w2 ~ treat_pmclong, final_drops), #long vs. PMC
  lm(matched_attriter_w2 ~ treat_vclong, final_drops), #long vs. VC
  lm(matched_attriter_w2 ~ treat_shortlong_collapse, final_drops), #long vs. both short videos
  lm(matched_attriter_w2 ~ short_pmc, final_drops), #short pmc vs. empty control
  lm(matched_attriter_w2 ~ short_vc, final_drops), #short vc vs. empty control
  lm(matched_attriter_w2 ~ treat_short_collapse, final_drops), #short vc vs. empty control
  lm(matched_attriter_w2 ~ pmc_vs_vc, final_drops) #short vs. short 
)

table(final_drops$pmc_vs_vc)

models_w3 <- list(
  lm(matched_attriter_w3 ~ treat_long_collapse2, final_drops), #long vs. placebo
  lm(matched_attriter_w3 ~ treat_long_collapse1, final_drops), #long vs. placebo+control 
  lm(matched_attriter_w3 ~ treat_pmclong, final_drops), #long vs. PMC
  lm(matched_attriter_w3 ~ treat_vclong, final_drops), #long vs. VC
  lm(matched_attriter_w3 ~ treat_shortlong_collapse, final_drops), #long vs. both short videos
  lm(matched_attriter_w3 ~ short_pmc, final_drops), #short pmc vs. empty control
  lm(matched_attriter_w3 ~ short_vc, final_drops), #short vc vs. empty control
  lm(matched_attriter_w3 ~ treat_short_collapse, final_drops), #short vc vs. empty control
  lm(matched_attriter_w3 ~ pmc_vs_vc, final_drops) #short vs. short 
)

models_w2_full <- list(
  lm(rand_attriter ~ treat_long_collapse2, final_drops), #long vs. placebo
  lm(rand_attriter ~ treat_long_collapse1, final_drops), #long vs. placebo+control 
  lm(rand_attriter ~ treat_pmclong, final_drops), #long vs. PMC
  lm(rand_attriter ~ treat_vclong, final_drops), #long vs. VC
  lm(rand_attriter ~ treat_shortlong_collapse, final_drops), #long vs. both short videos
  lm(rand_attriter ~ short_pmc, final_drops), #short pmc vs. empty control
  lm(rand_attriter ~ short_vc, final_drops), #short vc vs. empty control
  lm(rand_attriter ~ treat_short_collapse, final_drops), #short vc vs. empty control
  lm(rand_attriter ~ pmc_vs_vc, final_drops) #short vs. short 
)

models_w3_full <- list(
  lm(rand_attriter_w3 ~ treat_long_collapse2, final_drops), #long vs. placebo
  lm(rand_attriter_w3 ~ treat_long_collapse1, final_drops), #long vs. placebo+control 
  lm(rand_attriter_w3 ~ treat_pmclong, final_drops), #long vs. PMC
  lm(rand_attriter_w3 ~ treat_vclong, final_drops), #long vs. VC
  lm(rand_attriter_w3 ~ treat_shortlong_collapse, final_drops), #long vs. both short videos
  lm(rand_attriter_w3 ~ short_pmc, final_drops), #short pmc vs. empty control
  lm(rand_attriter_w3 ~ short_vc, final_drops), #short vc vs. empty control
  lm(rand_attriter_w3 ~ treat_short_collapse, final_drops), #short vc vs. empty control
  lm(rand_attriter_w3 ~ pmc_vs_vc, final_drops) #short vs. short 
)

# Create summary table
summary_table <- data.frame(
  Conditions = conditions,
  Wave2 = sapply(models_w2, extract_ols_results),
  Wave3 = sapply(models_w3, extract_ols_results),
  Wave2_FullSample = sapply(models_w2_full, extract_ols_results),
  Wave3_FullSample = sapply(models_w3_full, extract_ols_results)
)


means <- final_drops %>% 
  dplyr::group_by(video_treatment) %>% 
  summarize(
    n = n(),
    mean_w2 =mean(matched_attriter_w2),
    mean_w3 =mean(matched_attriter_w3),
    mean_w2_full = mean(rand_attriter),
    mean_w3_full = mean(rand_attriter_w3)
  )


# Print the table
print(summary_table)

# Create LaTeX table with xtable
latex_table <- xtable(
  summary_table, 
  caption = "Summary of Attrition Across Conditions and Waves",
  label = "tab:attrition_summary",
  align = c("l", "l", "c", "c", "c", "c")
)

# Specify column names
colnames(summary_table) <- c("Conditions", "Wave 2 (Matched)", "Wave 3 (Matched)", "Wave 2 (Full Sample)", "Wave 3 (Full Sample)")


# Pare down to five rows (Fill film vs. placebo, full film vs. placebo+control, full film vs. both short videos; Both Short vs. empty control; Short VC s. short PMC)

summary_table <- summary_table %>% 
  dplyr::filter(Conditions %in% c(
    "Full Film vs. Placebo",
    "Full Film vs. Placebo + Control",
    "Full Film vs. Both Short Videos",
    "Both Short Videos vs. Empty Control",
    "Short VC vs. Short PMC"
  ))

# Print LaTeX table to file
output_file <- "Figures/SI/attrition_summary.tex"

print(
  latex_table, 
  file = output_file, 
  include.rownames = FALSE,
  floating = TRUE,
  sanitize.text.function = identity, # To keep formatting like stars
  booktabs = TRUE, # For cleaner table style
  hline.after = c(-1, 0, nrow(summary_table)) # Horizontal lines at the top, header, and bottom
)

cat("LaTeX file saved as:", output_file, "\n")

output_file


#### Table SI.10 ####

#For this table, we manually moved the row with the count to the top and added an explanation column. We also renamed the variables so they are easier to interpret. While the formatting is a little different in the SI, the values are identical. 

#factor variables if needed...
final_drops <- final_drops %>% 
  mutate(
    region_2 = as.factor(region_2),
    vote_2016 = as.factor(vote_2016),
    vote_2020 = as.factor(vote_2020)
  )


# Convert factor variables to dummy variables and calculate means (if necessary)
final_drops_with_dummies <- final_drops %>%
  mutate(across(c(vote_2020, vote_2016, region_2), as.factor)) %>%
  clean_names() %>% # To ensure column names are compatible
  dummy_cols(select_columns = c("vote_2020", "vote_2016", "region_2"), remove_first_dummy = F) 

demos_with_dummies <- c("party_id", "ideology_yg", "age", "sex", "educ_dum", "white", "christian", "child", "job", "marr_dum", "turnout2020",
                        colnames(final_drops_with_dummies %>% select(starts_with("vote_2020_"), starts_with("vote_2016_"), starts_with("region_2_"))))


demos <- c("party_id", "ideology_yg", "age", "sex", "educ_dum", "white", "christian", "child", "job", "marr_dum", "turnout2020")


# Function to calculate mean, difference, SMD, and counts
calculate_smd <- function(group_1, group_2) {
  # Calculate mean and standard deviation for each group
  mean_1 <- mean(group_1, na.rm = TRUE)
  mean_2 <- mean(group_2, na.rm = TRUE)
  sd_1 <- sd(group_1, na.rm = TRUE)
  sd_2 <- sd(group_2, na.rm = TRUE)
  
  # Calculate difference in means
  difference <- mean_1 - mean_2
  
  # Calculate standardized mean difference (SMD)
  smd <- difference / sqrt((sd_1^2 + sd_2^2) / 2)
  
  # Calculate counts
  count_1 <- sum(!is.na(group_1))
  count_2 <- sum(!is.na(group_2))
  
  # Return the results as a list
  list(
    mean_1 = mean_1,
    mean_2 = mean_2,
    sd_1 = sd_1,
    sd_2 = sd_2,
    difference = difference,
    smd = smd,
    count_1 = count_1,
    count_2 = count_2
  )
}

# Filter the data to exclude NA values in treat_long_collapse2
filtered_data <- final_drops_with_dummies %>%
  filter(!is.na(treat_long_collapse2))



# Apply the function to every variable in demos_with_dummies
results <- map_dfr(demos, function(var) {
  group_1 <- filtered_data %>% filter(matched_attriter_w2 == FALSE) %>% pull(!!sym(var))
  group_2 <- filtered_data %>% filter(matched_attriter_w2 == TRUE) %>% pull(!!sym(var))
  
  smd_result <- calculate_smd(group_1, group_2)
  
  tibble(
    Variable = var,
    Mean_Non_Attriter = smd_result$mean_1,
    Mean_Attriter = smd_result$mean_2,
    Difference = smd_result$difference,
    SMD = smd_result$smd
  )
})

# Manually add the count row
count_row <- tibble(
  Variable = "Count",
  Mean_Non_Attriter = sum(!is.na(filtered_data %>% filter(matched_attriter_w2 == FALSE) %>% pull(demos_with_dummies[1]))),
  Mean_Attriter = sum(!is.na(filtered_data %>% filter(matched_attriter_w2 == TRUE) %>% pull(demos_with_dummies[1]))),
  Difference = NA,
  SMD = NA
)

# Bind the count row to the results
results <- bind_rows(results, count_row)

# Create LaTeX table
latex_table <- results %>%
  kbl(caption = "Comparison of Means for Attriters vs. Non-Attriters", format = "latex", booktabs = TRUE, digits = 2) %>%
  kable_styling(latex_options = c("hold_position"))

print(latex_table)


#### Table SI.11 ####

#create df for long vs. placebo (nationally representative)
long_placebo <- final_drops %>% 
  filter(!is.na(treat_long_collapse2))

#treatment assignment (w2 attrition)
attrition1 <- lm(matched_attriter_w2~treat_long_collapse2, long_placebo)


covariates_noimp <- paste("treat_long_collapse2", "partyID", "ideology_yg", "age", "sex", "educ_dum", "white", "region_2", "christian", "child", "job", "marr_dum", "turnout2020", sep = "+")

test1_noimp <- paste("matched_attriter_w2 ~ ", covariates_noimp, sep = "")

#model with covariates 
attrition1_covariates_noimp <- lm(test1_noimp, long_placebo)

####Interaction Terms####
interactions_noimp <-  paste("treat_long_collapse2*partyID", "ideology_yg*treat_long_collapse2", "age*treat_long_collapse2", "sex*treat_long_collapse2", "educ_dum*treat_long_collapse2", "white*treat_long_collapse2", "region_2*treat_long_collapse2", "christian*treat_long_collapse2", "child*treat_long_collapse2", "job*treat_long_collapse2", "marr_dum*treat_long_collapse2", "turnout2020*treat_long_collapse2", sep = "+")

test1_int_noimp <- paste("matched_attriter_w2 ~ ", interactions_noimp, sep = "")

#model with covariates + interactions
attrition1_interactions_noimp <- lm(test1_int_noimp, long_placebo)

attrition1_tables_noimp <- list(attrition1, attrition1_covariates_noimp, attrition1_interactions_noimp)

stargazer::stargazer(attrition1_tables_noimp, type="latex", title = "Likelihood of Attrition in Nationally Representative Sample: Non-Imputed Values (Long Video vs. Placebo)", dep.var.labels = "(Long vs. Placebo)", dep.var.caption = "Attriter", single.row = T)

```


## Main outcomes and analyses

```{r tables SI.12-15}

####Nationally Representative (SI.12)####
C1.1 <- lm_robust(aff_pol_idx~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

C1.2 <- lm_robust(aff_pol_idx_outparty~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

C1.3 <- lm_robust(BA_newsletter_clicked~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

C1.4 <- lm_robust(donate_any~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")


get_table_4(C1.1, C1.2, C1.3, C1.4, "Main Analysis: Nationally Representative Sample (Long vs. Placebo)", "Affective Polarization", "Outparty Only", "BA Newsletter", "Any Donation")

####no IPWs (SI.13)####
C1.1_noipws <- lm_robust(aff_pol_idx~treat_long_collapse2, final_drops, weights = weight_w2, se_type = "HC1")

C1.2_noipws <- lm_robust(aff_pol_idx_outparty~treat_long_collapse2, final_drops, weights = weight_w2, se_type = "HC1")

C1.3_noipws <- lm_robust(BA_newsletter_clicked~treat_long_collapse2, final_drops, weights = weight_w2, se_type = "HC1")

C1.4_noipws <- lm_robust(donate_any~treat_long_collapse2, final_drops, weights = weight_w2, se_type = "HC1")


get_table_4(C1.1_noipws, C1.2_noipws, C1.3_noipws, C1.4_noipws, "Main Analysis: Nationally Representative Sample (Long vs. Placebo; NO IPWs)", "Affective Polarization", "Outparty Only", "BA Newsletter", "Any Donation")


####Full Sample (SI.14)####
D1.1 <- lm_robust(aff_pol_idx~treat_long_collapse2, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

D1.2 <- lm_robust(aff_pol_idx_outparty~treat_long_collapse2, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

D1.3 <- lm_robust(BA_newsletter_clicked~treat_long_collapse2, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

D1.4 <- lm_robust(donate_any~treat_long_collapse2, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")


get_table_4(D1.1, D1.2, D1.3, D1.4, "Main Analysis: Full Sample (Long vs. Placebo)", "Affective Polarization", "Outparty Only", "BA Newsletter", "Any Donation")


####Full Sample no IPWs (SI.15)####
D1.1_noipw <- lm_robust(aff_pol_idx~treat_long_collapse2, final_drops, se_type = "HC1")

D1.2_noipw <- lm_robust(aff_pol_idx_outparty~treat_long_collapse2, final_drops, se_type = "HC1")

D1.3_noipw <- lm_robust(BA_newsletter_clicked~treat_long_collapse2, final_drops, se_type = "HC1")

D1.4_noipw <- lm_robust(donate_any~treat_long_collapse2, final_drops, se_type = "HC1")


get_table_4(D1.1_noipw, D1.2_noipw, D1.3_noipw, D1.4_noipw, "Main Analysis: Full Sample (Long vs. Placebo; NO IPWs)", "Affective Polarization", "Outparty Only", "BA Newsletter", "Any Donation")


```


## Donations extra analysis 

```{r tables SI.16-SI.18}

#create binary donation variables
final_drops <- final_drops %>% 
  mutate(
    donate_BA_any = if_else(donate_BA > 0, 1,0),
    donate_allsides_any = if_else(donate_allsides > 0, 1,0),
    donate_lrc_any = if_else(donate_lrc > 0, 1,0)
  )

####All Donations (Table SI.16)####
H1.2b_BA_t <- lm_robust(donate_BA~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

H1.2b_AS_t <- lm_robust(donate_allsides~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

H1.2b_LRC_t <- lm_robust(donate_lrc~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

H1.2b_all_t <- lm_robust(donate_all~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

get_table_4(H1.2b_BA_t, H1.2b_AS_t, H1.2b_LRC_t, H1.2b_all_t, "All Donations by Depolarization Organization", "Braver Angels", "AllSides", "LivingRoom Conversations", "Total")


#### NO IPWs (Table SI.17)####
H1.2b_BA_t_noIPWS <- lm_robust(donate_BA~treat_long_collapse2, final_drops, weights = weight_w2, se_type = "HC1")

H1.2b_AS_t_noIPWS <- lm_robust(donate_allsides~treat_long_collapse2, final_drops, weights = weight_w2, se_type = "HC1")

H1.2b_LRC_t_noIPWS <- lm_robust(donate_lrc~treat_long_collapse2, final_drops, weights = weight_w2, se_type = "HC1")

H1.2b_all_t_noIPWS <- lm_robust(donate_all~treat_long_collapse2, final_drops, weights = weight_w2, se_type = "HC1")

get_table_4(H1.2b_BA_t_noIPWS, H1.2b_AS_t_noIPWS, H1.2b_LRC_t_noIPWS, H1.2b_all_t_noIPWS, "All Donations by Depolarization Organization (NO IPWs)", "Braver Angels", "AllSides", "LivingRoom Conversations", "Total")


#### Any Donation (Table SI.18)

H1.2b_BA_any <- lm_robust(donate_BA_any~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

H1.2b_AS_any <- lm_robust(donate_allsides_any~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

H1.2b_LRC_any <- lm_robust(donate_lrc_any~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

get_table_3(H1.2b_BA_any, H1.2b_AS_any, H1.2b_LRC_any, "Any Donations by Depolarization Organization", "Braver Angels", "AllSides", "Living Room Conversations")


```


## Secondary Outcomes and Mechanisms

```{r tables SI.19-21}

####Optimism About Survival of Institutions: Table SI.19#### 
K1.1 <- lm_robust(optimism_survive_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

K1.2 <- lm_robust(optimism_survive_scale~treat_long_collapse2, final_drops, weights = weight_w2, se_type = "HC1")

K1.3 <- lm_robust(optimism_survive_scale~treat_long_collapse2, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

K1.4 <- lm_robust(optimism_survive_scale~treat_long_collapse2, final_drops, se_type = "HC1")

get_table_4(K1.1, K1.2, K1.3, K1.4, "Main Analysis: Optimism about Survival of Democratic Institutions Wave 2 (Long vs. Placebo)", "Matched (IPWs)", "Matched (No IPWs)", "Full (IPWs)", "Full (No IPWs)")


####Optimism about Restoring Civility (Table SI.20) ###
M1.1 <- lm_robust(optimism_civil_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

M1.2 <- lm_robust(optimism_civil_scale~treat_long_collapse2, final_drops, weights = weight_w2, se_type = "HC1")

M1.3 <- lm_robust(optimism_civil_scale~treat_long_collapse2, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

M1.4 <- lm_robust(optimism_civil_scale~treat_long_collapse2, final_drops, se_type = "HC1")

get_table_4(M1.1, M1.2, M1.3, M1.4, "Main Analysis: Optimism about Restoring Civility Wave 2 (Long vs. Placebo)", "Matched (IPWs)", "Matched (No IPWs)", "Full (IPWs)", "Full (No IPWs)")

####Dialogue as effective tool for change (SI.21) ####
N1.1 <- lm_robust(dialogue_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

N1.2 <- lm_robust(dialogue_scale~treat_long_collapse2, final_drops, weights = weight_w2, se_type = "HC1")

N1.3 <- lm_robust(dialogue_scale~treat_long_collapse2, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

N1.4 <- lm_robust(dialogue_scale~treat_long_collapse2, final_drops, se_type = "HC1")

get_table_4(N1.1, N1.2, N1.3, N1.4, "Main Analysis: Dialogue as Effective Tool For Change Wave 2 (Long vs. Placebo)", "Matched (IPWs)", "Matched (No IPWs)", "Full (IPWs)", "Full (No IPWs)")


#### Anti-Democratic Attitudes (SI.22) ####

#ipws
E1.1 <- lm_robust(antidem_idx~treat_long_collapse2, final_drops, se_type = "HC1", weights = ipw_matched_w2_noimp)

#no IPWs
E1.2 <- lm_robust(antidem_idx~treat_long_collapse2, final_drops, se_type = "HC1", weights = weight_w2)

#full sample
E1.3 <- lm_robust(antidem_idx~treat_long_collapse2, final_drops, se_type = "HC1", weights = ipw_full_w2_noimp)

#no IPWs
E1.4 <- lm_robust(antidem_idx~treat_long_collapse2, final_drops, se_type = "HC1")

get_table_4(E1.1, E1.2, E1.3,E1.4, "Main Analysis: Anti-Democratic Attitudes Wave 2 (Long vs. Placebo)", "Matched Sample (IPWs)", "Matched Sample (No IPWs)", "Full Sample (IPWs)", "Full Sample (No IPWs")


```

```{r Possible mechanisms Figure SI.8}

#higher more united
final_drops <- final_drops %>% 
  mutate(
    division = as.numeric(Q52_w2),
    division= if_else(division>100, NA, division),
    division_scaled = scale(division)
  )

#Regression models 
med_mass.perc <- lm(mass_perc_idx~treat_long_collapse2, final_drops, weights=ipw_matched_w2_noimp)

med_division <- lm(division_scaled~treat_long_collapse2, final_drops, weights=ipw_matched_w2_noimp)

med_openmind <- lm(understand_total~treat_long_collapse2, final_drops, weights=ipw_matched_w2_noimp)

med_empathy <- lm(empathy_scale~treat_long_collapse2, final_drops, weights=ipw_matched_w2_noimp)

stereo_neg <- lm(stereo_outparty_neg_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)

stereo_pos <- lm(stereo_outparty_pos_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)


#Plot
modelslist <- list(med_empathy, med_openmind, med_division, med_mass.perc, stereo_neg, stereo_pos)

nameslist <- as.character(c("Empathy for outgroup", "Open-mindedness", "Perceived division", "Mass perceptions", "Negative stereotypes", "Positive stereotypes"))

mediators <- plot_generate_6(med_empathy, med_openmind, med_division, med_mass.perc, stereo_neg, stereo_pos, "Empathy for outgroup", "Open-mindedness", "Perceived unity", "Mass perceptions", "Negative stereotypes", "Positive stereotypes", "", c(-0.5,0.5))

mediators

#save
ggsave("Figures/SI/Figure_SI.8.png", plot = mediators, width=10, height=7)


```

## Extensions 

```{r HTEs (Tables SI.23-25)}

## Note that treat_long_collapse2 is the binary treatment indicator. In the manuscript, the name is manually changed to Full Film for clarity of presentation. 

#####Party ID (Table SI.23)####
O1.1 <- lm_robust(aff_pol_idx~treat_long_collapse2*partyID, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

O1.2 <- lm_robust(aff_pol_idx~treat_long_collapse2*partyID, final_drops, weights = weight_w2, se_type = "HC1")

O1.3 <- lm_robust(aff_pol_idx~treat_long_collapse2*partyID, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

O1.4 <- lm_robust(aff_pol_idx~treat_long_collapse2*partyID, final_drops, se_type = "HC1")

get_table_4(O1.1, O1.2, O1.3, O1.4, "Main Analysis: Affective Polarization HTE by Party ID Wave 2 (Long vs. Placebo)", "Matched (IPWs)", "Matched (No IPWs)", "Full (IPWs)", "Full (No IPWs)")


####Ideology (Table SI.24)####
P1.1 <- lm_robust(aff_pol_idx~treat_long_collapse2*ideo, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

P1.2 <- lm_robust(aff_pol_idx~treat_long_collapse2*ideo, final_drops, weights = weight_w2, se_type = "HC1")

P1.3 <- lm_robust(aff_pol_idx~treat_long_collapse2*ideo, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

P1.4 <- lm_robust(aff_pol_idx~treat_long_collapse2*ideo, final_drops, se_type = "HC1")

get_table_4(P1.1, P1.2, P1.3, P1.4, "Main Analysis: Affective Polarization HTE by Ideology Wave 2 (Long vs. Placebo)", "Matched (IPWs)", "Matched (No IPWs)", "Full (IPWs)", "Full (No IPWs)")


#### Media Index (Table SI.26)####
# media_truth_dum is manually changed to Media Trust in the SI for ease of interpretation 

Q1.1 <- lm_robust(aff_pol_idx~treat_long_collapse2*media_truth_dum, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

Q1.2 <- lm_robust(aff_pol_idx~treat_long_collapse2*media_truth_dum, final_drops, weights = weight_w2, se_type = "HC1")

Q1.3 <- lm_robust(aff_pol_idx~treat_long_collapse2*media_truth_dum, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

Q1.4 <- lm_robust(aff_pol_idx~treat_long_collapse2*media_truth_dum, final_drops, se_type = "HC1")

get_table_4(Q1.1, Q1.2, Q1.3, Q1.4, "Main Analysis: Affective Polarization HTE by Partisan Media Confidence Wave 2 (Long vs. Placebo)", "Matched (IPWs)", "Matched (No IPWs)", "Full (IPWs)", "Full (No IPWs)")


```

```{r wave 3 (Tables SI.26-27)}

#Table SI.26
R1.1 <- lm_robust(aff_pol_idx_w3~treat_long_collapse2, final_drops, weights = ipw_matched_w3_noimp, se_type = "HC1")

R1.2 <- lm_robust(aff_pol_idx_w3~treat_long_collapse2, final_drops, weights = weight_w3, se_type = "HC1")

R1.3 <- lm_robust(aff_pol_idx_w3~treat_long_collapse2, final_drops, weights = ipw_full_w3_noimp, se_type = "HC1")

R1.4 <- lm_robust(aff_pol_idx_w3~treat_long_collapse2, final_drops, se_type = "HC1")

get_table_4(R1.1, R1.2, R1.3,R1.4, "Main Analysis: Affective Polarization Wave 3 (Long vs. Placebo)", "Matched (IPWs)", "Matched (No IPWs)", "Full (IPWs)", "Full (No IPWs)")

#Table SI.26
S1.1 <- lm_robust(aff_pol_idx_outparty_w3~treat_long_collapse2, final_drops, weights = ipw_matched_w3_noimp, se_type = "HC1")

S1.2 <- lm_robust(aff_pol_idx_outparty_w3~treat_long_collapse2, final_drops, weights = weight_w3, se_type = "HC1")

S1.3 <- lm_robust(aff_pol_idx_outparty_w3~treat_long_collapse2, final_drops, weights = ipw_full_w3_noimp, se_type = "HC1")

S1.4 <- lm_robust(aff_pol_idx_outparty_w3~treat_long_collapse2, final_drops, se_type = "HC1")

get_table_4(S1.1, S1.2, S1.3,S1.4, "Main Analysis: Affective Polarization Outparty Only Wave 3 (Long vs. Placebo)", "Matched (IPWs)", "Matched (No IPWs)", "Full (IPWs)", "Full (No IPWs)")


```

## Pre-registered analyses not featured in manuscript

```{r tables SI.28-34}

####Table SI.28####
A1.1 <- lm_robust(aff_pol_idx~treat_long_collapse1, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

summary(A1.1)

A1.2 <- lm_robust(aff_pol_idx_outparty~treat_long_collapse1, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

summary(A1.2)

A1.3 <- lm_robust(BA_newsletter_clicked~treat_long_collapse1, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

summary(A1.3)

A1.4 <- lm_robust(donate_any~treat_long_collapse1, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

summary(A1.4)

get_table_4(A1.1, A1.2, A1.3, A1.4, "Main Analysis: Nationally Representative Sample Wave 2 (Long vs. Placebo+Control)", "Affective Polarization", "Outparty Only", "BA Newsletter", "Any Donation")

####Table SI.29####
B1.1 <- lm_robust(aff_pol_idx~treat_long_collapse1, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")


B1.2 <- lm_robust(aff_pol_idx_outparty~treat_long_collapse1, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")


B1.3 <- lm_robust(BA_newsletter_clicked~treat_long_collapse1, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")


B1.4 <- lm_robust(donate_any~treat_long_collapse1, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

get_table_4(B1.1, B1.2, B1.3, B1.4, "Main Analysis: Full Sample Wave 2 (Long vs. Placebo+Control)", "Affective Polarization", "Outparty Only", "BA Newsletter", "Any Donation")


####Table SI.30####
#ensure factor varialbe
final_drops <- final_drops %>% 
  mutate(
    short_vs_long = as.factor(treat_long_short)
  )

#0 is empty control; 1 is Short PMC; 2 is Short VC; 3 is full film. 

#Short vs. long 1
U1.1 <- lm_robust(aff_pol_idx~short_vs_long, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")


U1.2 <- lm_robust(aff_pol_idx_outparty~short_vs_long, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")


U1.3 <- lm_robust(BA_newsletter_clicked~short_vs_long, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")


U1.4 <- lm_robust(donate_any~short_vs_long, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

get_table_4(U1.1, U1.2, U1.3, U1.4, "Main Analysis: Nationally Representative Sample  Wave 2 (Full Film vs. Short Videos)", "Affective Polarization", "Outparty Only", "BA Newsletter", "Any Donation")


#### Table SI.31####

T1.1 <- lm_robust(aff_pol_idx~short_vs_long, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")


T1.2 <- lm_robust(aff_pol_idx_outparty~short_vs_long, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")


T1.3 <- lm_robust(BA_newsletter_clicked~short_vs_long, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")


T1.4 <- lm_robust(donate_any~short_vs_long, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

get_table_4(T1.1, T1.2, T1.3, T1.4, "Main Analysis: Full Sample  Wave 2 (Full Film vs. Short Videos)", "Affective Polarization", "Outparty Only", "BA Newsletter", "Any Donation")


#### Table SI.32 ####
R2.1 <- lm_robust(aff_pol_idx~treat_short_collapse, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

R2.2 <- lm_robust(aff_pol_idx_outparty~treat_short_collapse, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

R2.3 <- lm_robust(BA_newsletter_clicked~treat_short_collapse, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

R2.4 <- lm_robust(donate_any~treat_short_collapse, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

get_table_4(R2.1, R2.2, R2.3, R2.4, "Excluded Pre-Specified Analysis: Short Films Main Outcomes Wave 2 (Pooled short vs. Empty control)", "Affective Polarization", "Outparty Affective Polarization", "BA Newsletter Clicks", "Any Donation")

#### Table SI.33 #### 
V1.1 <- lm_robust(UA_interest~treat_long_collapse2, final_drops, weights = ipw_matched_w3_noimp, se_type = "HC1")

V1.2 <- lm_robust(UA_interest~treat_long_collapse2, final_drops, weights = weight_w3, se_type = "HC1")

V1.3 <- lm_robust(UA_interest~treat_long_collapse2, final_drops, weights = ipw_full_w3_noimp, se_type = "HC1")

V1.4 <- lm_robust(UA_interest~treat_long_collapse2, final_drops, se_type = "HC1")

get_table_4(V1.1, V1.2, V1.3, V1.4, "Excluded Pre-Specified Analysis: Unify America Sign-Ups (Full Film vs. Placebo)", "Matched (IPWs)", "Matched (No IPWs)", "Full (IPWs)", "Full (No IPWs)")

#### Table SI.34 #### 
W1.1 <- lm_robust(antidem_idx_w3~treat_long_collapse2, final_drops, weights=ipw_matched_w3_noimp, se_type = "HC1")

W1.2 <- lm_robust(antidem_idx_w3~treat_long_collapse2, final_drops, weights=weight_w3, se_type = "HC1")

W1.3 <- lm_robust(antidem_idx_w3~treat_long_collapse2, final_drops, weights=ipw_full_w3_noimp, se_type = "HC1")

W1.4 <- lm_robust(antidem_idx_w3~treat_long_collapse2, final_drops, se_type = "HC1")

get_table_4(W1.1, W1.2, W1.3, W1.4, "Excluded Pre-Specified Analysis: Anti-Democratic Attitudes Wave 3 (Full Film vs. Placebo)", "Matched (IPWs)", "Matched (No IPWs)", "Full (IPWs)", "Full (No IPWs)")


```

## Ancillary Analyses

```{r Affective Polarization Index Breakdown}

#### Figure SI.9 ####
therm <- lm(therm_diff_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)

trust <- lm(trust_diff_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)

discomfort_marr <- lm(discomfort_outparty_marriage_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)

discomfort_friend <- lm(discomfort_outparty_friends_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)

discomfort_neighbor <- lm(discomfort_outparty_neighbors_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)

threat <- lm(threat_outparty_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)

neg_par <- lm(neg_partisan_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)


models_list <- list(therm, trust, discomfort_marr, discomfort_friend, discomfort_neighbor, threat, neg_par)

names_list <- c("Feeling Thermometer", "Trust Difference",
               "Discomfort Marriage", "Discomfort
               Friends", "Discomfort Neighbors", 
               "Outparty Threat", "Negative Partisanship")

affpol_components <- plot_generate_7b(models_list, names_list, "Affective Polarization Index Breakdown by Component", c(-0.6, 0.4)) + theme_minimal()

affpol_components

ggsave("Figures/SI/Figure_SI.9.png", plot = affpol_components, width=10, height=7)


#### Figure SI.10 #### 
#Outparty Only
therm_out <- lm(therm_outparty_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)

trust_out <- lm(trust_outparty_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp)


models_list2 <- list(therm_out, trust_out, discomfort_marr, discomfort_friend, discomfort_neighbor, threat, neg_par)

names_list2 <- c("Feeling Thermometer Outparty", "Trust Outparty","Discomfort Marriage", "Discomfort Friends", "Discomfort Neighbors", "Outparty Threat", "Negative Partisanship")


affpol_components_outparty <- plot_generate_7b(models_list2, names_list2, "Affective Polarization Index (Outparty Only): Breakdown by Component", c(-0.6, 0.4)) + theme_minimal()

affpol_components_outparty

ggsave("Figures/SI/Figure_SI.10.png", plot = affpol_components_outparty, width=10, height=7)


#### Table SI.35 #### 
final_drops <- final_drops %>% 
  mutate(
    therm_inparty_scale = scale(therm_inparty),
    therm_inparty_w3_scale = scale(therm_inparty_w3),
    trust_inparty_scale = scale(trust_inparty),
    trust_inparty_w3_scale = scale(trust_inparty_w3)
  )

trust_inparty_only <- lm_robust(trust_inparty_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

therm_inparty_only <- lm_robust(therm_inparty_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")


get_table_2(trust_inparty_only, therm_inparty_only, "Effect of Documentary on Ingroup Affective Polarization Measures","Ingroup Trust", "Ingroup Warmth")


```


```{r HTEs secondary outcomes}

#### Figure SI.11 #### 
HTE1.3_dialogue <- lm(dialogue_scale~treat_long_collapse2*partyID, final_drops, weights = ipw_matched_w2_noimp)

#generate df and plot
HTE1.3_df <- get_cate(HTE1.3_dialogue, "Democrats", "Republicans")
dialogue_cate <- plot_cate(HTE1.3_df, "A: Dialogue Effectiveness", c(-0.5,1))

dialogue_cate

#optimism about restoring civility
HTE1.5_optimism_civil <- lm(optimism_civil_scale~treat_long_collapse2*partyID, final_drops, weights = ipw_matched_w2_noimp)

#generate df and plot
HTE1.5_df <- get_cate(HTE1.5_optimism_civil, "Democrats", "Republicans")
optimism_civil_cate <- plot_cate(HTE1.5_df, "B: Restoring Civility and Goodwill", c(-0.5,1))

optimism_civil_cate

#optimism about survival of democratic institutions 
HTE1.6_optimism_survive <- lm(optimism_survive_scale~treat_long_collapse2*partyID, final_drops, weights = ipw_matched_w2_noimp)

#generate df and plot
HTE1.6_df <- get_cate(HTE1.6_optimism_survive, "Democrats", "Republicans")
optimism_democracy_cate <- plot_cate(HTE1.6_df, "C: Survival of Democratic Institutions", c(-0.5,1))

optimism_democracy_cate

#non-violent change still possible 
HTE1.8_nonviol <- lm(nonviol_change_scale~treat_long_collapse2*partyID, final_drops, weights = ipw_matched_w2_noimp)

#generate df and plot
HTE1.8_df <- get_cate(HTE1.8_nonviol, "Democrats", "Republicans")
nonviol_cate <- plot_cate(HTE1.8_df, "D: Non-Violent Change", c(-1,1))

nonviol_cate

library(patchwork)

combined_secondary_HTE_plots <- dialogue_cate + optimism_civil_cate+ optimism_democracy_cate + nonviol_cate + plot_layout(guides = "collect")

combined_secondary_HTE_plots

#save
ggsave("Figures/SI/Figure_SI.12.png", plot = combined_secondary_HTE_plots, width=10, height=7)


#### Figure SI.12 #### 
antidem_cate_ID <- lm(antidem_idx~treat_long_collapse2*partyID, final_drops, weights = ipw_matched_w2_noimp)

antidem_partyID <- get_cate(antidem_cate_ID, "Democrats", "Republicans")

antidem_partyID_plot <- plot_cate(antidem_partyID, "CATE Anti-Democratic Attitudes Party ID:
          Nationally Representative Sample Wave 2", c(-0.5,0.5))

antidem_partyID_plot

ggsave("Figures/SI/Figure_SI.12.png", plot = antidem_partyID_plot, width=10, height=7)


```

```{r survey response time}

#### Table SI.36 ####

# Create time bins by 24 hours
final_drops <- final_drops %>% 
  mutate(
    time_bins = case_when(
      post_treat_w2 <= 24 ~ 0,  # First bin: 24 hours or less
      post_treat_w2 > 24 & post_treat_w2 <= 48 ~ 1,  # Second bin: 25-48 hours
      post_treat_w2 > 48 & post_treat_w2 <= 72 ~ 2,  # Third bin: 49-72 hours
      post_treat_w2 > 72 & post_treat_w2 <= 96 ~ 3,  # Fourth bin: 73-96 hours
      post_treat_w2 > 96 & post_treat_w2 <= 120 ~ 4, # Fifth bin: 97-120 hours
      post_treat_w2 > 120 & post_treat_w2 <= 144 ~ 5, # Sixth bin: 121-144 hours
      post_treat_w2 > 144 & post_treat_w2 <= 168 ~ 6, # Seventh bin: 145-168 hours
      post_treat_w2 > 168 & post_treat_w2 <= 192 ~ 7, # Eighth bin: 169-192 hours
      post_treat_w2 > 192 & post_treat_w2 <= 216 ~ 8, # Ninth bin: 193-216 hours
      post_treat_w2 > 216 ~ 9  # Tenth bin: 217 hours or more
    )
  )

by_time <- lm_robust(aff_pol_idx~treat_long_collapse2*time_bins, final_drops, weights=ipw_matched_w2_noimp, se_type = "HC1")

summary(by_time)

by_time_out <- lm_robust(aff_pol_idx_outparty~treat_long_collapse2*time_bins, final_drops, weights=ipw_matched_w2_noimp, se_type = "HC1")

summary(by_time_out)

#make table
get_table_2(by_time, by_time_out, "Affective Polarization by Wave 2 Response Time", "Affective Polarization", "Outparty Only")

```


```{r longitudinal secondary outcomes}

#### Figure SI.13 ####
#belief in possibility of non-violent democratic change 
Y1.1_w3 <- lm(nonviol_change_scale_w3~treat_long_collapse2, final_drops, weights = ipw_matched_w3_noimp)

Y1.4_w3 <- lm(dialogue_scale_w3~treat_long_collapse2, final_drops, weights = ipw_matched_w3_noimp)

Y2.1_w3 <- lm(optimism_civil_scale_w3~treat_long_collapse2, final_drops, weights = ipw_matched_w3_noimp)

Y2.2_w3 <- lm(optimism_survive_scale_w3~treat_long_collapse2, final_drops, weights = ipw_matched_w3_noimp)


optimism_plot_PAP_w3 <- plot_generate_4(Y1.1_w3, Y1.4_w3, Y2.1_w3, Y2.2_w3, "Non-Violent Change Still Possible", "Belief in Efficacy of Dialogue","Restoring Civility", "Survival of Democratic Institutions","", c(-0.3, .6))

print(optimism_plot_PAP_w3)

#save
ggsave("Figures/SI/Figure_SI.13.png", plot = optimism_plot_PAP_w3, width=10, height=7)


```

```{r HTE with covars}

#### Table SI.37 #### 
#HTEs covariates 
covars <- c("age", "sex", "educ_dum", "white", "christian", "child", "job", "marr_dum", "turnout2020", "ideo_scale")
covars_formula <- paste(covars, collapse = " + ")
formula <- as.formula(paste("aff_pol_idx ~ treat_long_collapse2 * partyID +", covars_formula))

HTE1.1_covars <- lm_robust(formula, data = final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

HTE1.1_covars_2 <- lm(formula, data = final_drops, weights = ipw_matched_w2_noimp)

summary(HTE1.1_covars_2)

#HTEs covariates outparty 
formula_out <- as.formula(paste("aff_pol_idx_outparty ~ treat_long_collapse2 * partyID +", covars_formula))

HTE1.1_covars_out <- lm_robust(formula_out, data = final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")
summary(HTE1.1_covars_out)

HTE1.1_covars_out_2 <- lm(formula_out, data = final_drops, weights = ipw_matched_w2_noimp)


#create table 
get_table_2(HTE1.1_covars, HTE1.1_covars_out, "", "Affective polarization", "Outparty only")


```

```{r additional secondary outcomes not pre-registered}

#### Table SI.38 ####
L1.1 <- lm_robust(optimism_pol_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

L1.2 <- lm_robust(optimism_pol_scale~treat_long_collapse2, final_drops, weights = weight_w2, se_type = "HC1")

L1.3 <- lm_robust(optimism_pol_scale~treat_long_collapse2, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

L1.4 <- lm_robust(optimism_pol_scale~treat_long_collapse2, final_drops, se_type = "HC1")

get_table_4(L1.1, L1.2, L1.3, L1.4, "Main Analysis: Optimism about Overcoming Polarization Wave 2 (Long vs. Placebo)", "Matched (IPWs)", "Matched (No IPWs)", "Full (IPWs)", "Full (No IPWs)")

#### Table SI.39 #### 
Y1.1 <- lm_robust(rebuild_trust_scale~treat_long_collapse2, final_drops, weights = ipw_matched_w2_noimp, se_type = "HC1")

Y1.2 <- lm_robust(rebuild_trust_scale~treat_long_collapse2, final_drops, weights = weight_w2, se_type = "HC1")

Y1.3 <- lm_robust(rebuild_trust_scale~treat_long_collapse2, final_drops, weights = ipw_full_w2_noimp, se_type = "HC1")

Y1.4 <- lm_robust(rebuild_trust_scale~treat_long_collapse2, final_drops, se_type = "HC1")

get_table_4(Y1.1, Y1.2, Y1.3, Y1.4, "Main Analysis: Time It Would Take to Rebuild Trust (reverse-coded) Wave 2 (Long vs. Placebo)", "Matched (IPWs)", "Matched (No IPWs)", "Full (IPWs)", "Full (No IPWs)")


```

```{r comparison betwen short videos and full documentary}

#### Figure SI.14 ####
#both short videos vs. long video 
final_drops <- final_drops %>% 
  mutate(
    treat_shortlong_collapse =
      case_when(
      video_treatment == "Treatment Long" ~ 1,
      video_treatment == "Treatment Short PMC" ~ 0,
      video_treatment == "Treatment Short VC" ~ 0
    )
  )

#polarization 
full_vs_short <- lm(aff_pol_idx~treat_shortlong_collapse, final_drops, weights = ipw_matched_w2_noimp)

#outparty 
full_vs_short_out <- lm(aff_pol_idx_outparty~treat_shortlong_collapse, final_drops, weights = ipw_matched_w2_noimp)


#BA newsletter
full_vs_short_newsletter <- lm(BA_newsletter_clicked~treat_shortlong_collapse, final_drops, weights = ipw_matched_w2_noimp)

#donation
full_vs_short_donate <- lm(donate_any~treat_shortlong_collapse, final_drops, weights = ipw_matched_w2_noimp)


short_vs_long_plot <- plot_generate_4_2(full_vs_short, full_vs_short_out, full_vs_short_newsletter, full_vs_short_donate, "Affective Polarization", "Outparty Only", "BA Newsletter", "Any Donation", "", c(-.2, .2))

short_vs_long_plot

#save
ggsave("Figures/SI/Figure_SI.14.png", plot = short_vs_long_plot, width=10, height=7)

```

